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Unmanned Air Vehicles have become increasingly important on the modern 
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fidelity, non-linear equations of motion were derived in matrix form that represented 
any six degree of freedom aircraft model, and were then tailored for use on specific 
aircraft. Computer modeling of the resulting equations of motion, as well as the 
sensors used on the aircraft, was done using SIMULINK and MATLAB software. The 
resulting computer model provides a non-linear system of equations, which are easily 
linearized at any desired flight condition, as required by the proposed control and 
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I. INTRODUCTION 


A. IMPORTANCE OF UNMANNED AIR VEHICLES 


Unmanned Aerial Vehicles have become increasingly more important, both on 
the battlefield and in civilian service, since the Ryan Q-2C Firebee target drone 
introduced the “modern age” of UAVs in 1960 [Ref. Siu 91]. From that time on, 
military planners have assured that UAVs have the capability to collect intelligence, 
target enemy positions, gather bomb damage assessment, as well as perform many 
other tasks. The real benefit in using unmanned aircraft lies in the fact that many 
missions can be performed deep in enemy territory, all without endangering the lives 
of pilots, or risking the loss of a much more expensive aircraft. With the recent use of 
UAVs in Operation Desert Storm, improvements in the current technology are both 
indicated and desirable. 

The most capable UAV in service of the United States Navy and Marine Corps 
today is the Pioneer Short Range UAV. The system is hampered, though, by the 
large amount of runway and special equipment needed to launch and land the aircraft. 
These requirements limit the usefulness of tHe Rionéer by keeping the aircraft take- 
off and landing area well away from the areas where the ground forces are operating. 
This distance then leads to longer transit times to and from the assigned operating 
area, and thus a shorter time on station. However, what is really required in many 
instances by the ground forces is an aircraft. that can respond quickly to a changing 
tactical situation. The Airborne Remotely Operated Device, AROD, was an attempt 
by Sandia National Laboratories [Ref. Wh 87], in response to a requirement by the 


Naval Ocean Systems Center, NOSC, to respond to these needs. 


The United States Marine Corps had set a requirement for a short range, direct 


support UAV as described in [Ref. MCG 87]: 


“...to allow the front line commander to see “over the next hill”, to a 


distance of two kilometers...” 


The AROD was designed to be a ducted fan, hovering device carrying a fiber optic 
data link and on board cameras. AROD testing was canceled [Ref. Sa 89] as the 
Department of Defense requirements grew, requiring a minimum range of 30 km for 
all Short Range UAVs. The AROD was incapable of this kind of range, however; the 
design is still potentially useful. 

The Unmanned Air Vehicle Flight Research Lab, UAV FRL, at the Naval Post- 
graduate School has proposed a solution using the AROD that would satisfy the 
DoD short range UAV requirements, while maintaining the important capability for 


vertical take-offs and landings. 


B. UAV RESEARCH USING THE AROD VEHICLE 

Unmanned Aerial Vehicle research underway at NPS has taken the AROD air- 
frame and fitted it with wings from the Aquila UAV [Ref. Kre 92, Sto 93] in order 
to give the AROD forward flight capability. The proposed configuration will give the 
Archytas (an AROD with wings) the ability to take-off and land vertically and then 
transition to horizontal flight for the mission. This design will explore new technol- 
ogy, driven by the goals established by the UAV Joint Project Office [Ref. DOD 92] 
of: 


e Take off weight under 200 Ib 
e Carry a 50 lb payload 


e Fly at a maximum speed of 150 kts 


to 


e Take-off and land in an area 30m by 60m 


The vertical flight 1s accomplished with a powerful ducted fan,which causes a great 
deal of gyroscopic coupling and torque when producing enough thrust to lift the 
aircraft. Therefore, in order to achieve stable take—offs and landings, a three—axis 
autopilot is a necessary feature of the aircraft. Additional capabilities desired in 
the final version of the Archytas are guidance and navigation systems which will 
allow autonomous operation, as well as a Global Positioning System aided autoland 


capability. 


C. REQUIREMENT FOR MODELING 


Simulation and modeling of the aircraft are essential to the successful design 
of a control system capable of autonomous flight. The model must be a very high 
fidelity, non-linear model, that can be easily linearized at any given flight condition. 
The model should be able to interpolate between data points resulting from wind 
tunnel testing in order to simulate the highly non-linear transition from vertical to 
horizontal flight. Moreover, the model must also be capable of including the outputs 
of the sensors as inputs to the control and navigation system for sensors located at 
any arbitrary location on the aircraft. 

This thesis develops a six degree of freedom model for the AROD in the vertical 
flight regime, as well as for an aircraft in a fixed wing configuration. This test aircraft, 
named Bluebird, is used to test the guidance, navigation, and control, GNC, systems 
in horizontal flight, since there currently is no aerodynamic data available for the 
Archytas configuration. Use of the Bluebird will provide the capability to design and 


test a GNC system on a stable aircraft before the first flight on the Archytas. 


Il. BACKGROUND 


A. DESCRIPTION OF AROD 

The AROD was designed by the Sandia Research Laboratory in Albuquerque, 
New Mexico in a project managed by NOSC. The vehicle possessed no flying surfaces 
and relied solely on powered lift for flight. Control of the aircraft was obtained through 
the use of four fixed anti-torque vanes and four moveable control vanes positioned in 
the propeller wash of the duct [Ref. We 88]. The main features of the AROD were 
vertical take-off and landing, VTOL, flight, lightweight construction, compact size, 
and minimal support equipment required. However, the AROD required most of the 
engine output to maintain the powered lift, so very little excess thrust was left for 
translational flight. 

An important aspect of the AROD design was the improvement in static perfor- 
mance provided by the efficiency of the ducted fan design. The addition of the shroud 
around the three—bladed propeller resulted in increased mass flow through the fan, 
and thus more static thrust when compared to a conventional propeller configuration 
(Ref. Kre 92]. The AROD is shown in Figure 2.1 and characteristics of the AROD 
are tabulated in Table 2.1. The moveable control vanes are al] used in combination 
to exert the desired control forces on AROD. Roll control is achieved by deflecting the 
four vanes in the same direction, while pitch and yaw control is obtained by deflecting 
a pair of vanes in the required direction. The numbering of the vanes is shown in 
Figure 2.2 and the combinations of vane deflections required for positive roll, pitch, 


and yaw motion are given in Table 2.2. 





Figure 2.1: Airborne Remotely Operated Device, (Ref. Siu 91] 


B. AROD MODELING 

The AROD vehicle has been the subject of several theses at NPS. The designs 
presented in the theses rely on the AROD model given by Sandia Labs in the original 
design of their controller [Ref. Wh 87, Wh 91]. This model was based on the more 
classical technique of linearizing an aircraft model, based on dimensional derivatives 
in a state space form. The resulting model was acceptable for the AROD in a hovering 
and near vertical translational flight mode, but was not easily adaptable to anything 
other than the narrow range of conditions planned for AROD. The Sandia Labs 


papers also pointed out several types of coupling in the AROD. The most prominent 


TABLE 2.1: PHYSICAL CHARACTERISTICS OF AROD 


[thet Diameter, AY 
[Propeller Radius, Rin] 
[Ext Radus 37 
[——Tilet Area Ratio 1219 J] 
[Exit Area Ratio [1.15 
[Exterior Contour | Tapered Rear] 
[Propeller Location, % chord [25% | 
[Number of Blades | 3 J 
[Engine Speed, Max | 8000 pm] 
[Engine Speed, Nom. | 6500 rpm | 
[Tip Speed, Max. | 838 lpm] 
[Tp Speed, Nom. | 680-fpm | 
[Mass Moment of Inertia, Zz | 1231 slug 77] 
[Mass Moment of Inertia, I, | 3.9548 slug = 7] 
[Mass Moment of Inertia. I, | 3.9825 slug =| 
[ Prop Mass Moment of Inertia, Ir | 0.00898 slug = 7] 
] Prop Mass Moment of Inertia, J,, | 0.0045 slug — f? | 
[ Prop Mass Moment of Inertia. . | 0.0005 slug=f™ | 























TABLE 2.2: VANE DEFLECTION COMBINATIONS FOR POSITIVE 
ANGLES 


|__| Vane Combination | 
| Roll, @ | M+e+ +N | 


[PichO|  -™ | 
[Yaw vf hI 





of the coupling effects is the gyroscopic coupling between the pitch and yaw axes 
resulting from the large amount of angular momentum contributed to the aircraft 
by the propeller. Another dynamic coupling exists between the altitude-rate and the 
vehicle attitude, since a loss of lift due to thrust will occur when the vehicle is tilted to 


generate horizontal motion. Yet a third dynamic coupling exists between the altitude 
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Figure 2.2: AROD Direction of Positive Vane Deflections 


and roll control loops, since the reactive torques applied to the roll axis vary as the 
engine speed is varied. Sandia Labs also provided data for modeling both the engine 
and the servos as second order transfer functions which were used in this thesis. 
Additional information was obtained by Weir [Ref. We 88] in wind tunnel test- 
ing. This information included non-dimensional derivatives for vane effectiveness and 
non-dimensional stability derivatives. The report also stated that the control-vane 
effectiveness is constant out to at least 25 deg of deflection. Wind tunnel data were 
also presented to show that control vane effectiveness in approximately the same for 
translational flight as for hovering flight. This equivalence is due to the fact that the 
vanes are located in the high speed flow aft of the propeller and are not significantly 


affected by the freestream. 


TABLE 2.3: PHYSICAL CHARACTERISTICS OF BLUEBIRD 


[_Wegh —s«dYS~Ci*é“‘( TS 
[Average Wing Chord, @ | 1.802 f__| 
[Wing Span,b | aa f__| 
[ Planform Surface Area, S| 22.380 f7_| 
[Engine Power | 40 HP _| 
[ Mass Moment of Inertia, I, 100 slug —J* | 
[ Mass Moment of Inertia, Ty | 16.12 slug — J*| 
[ Mass Moment of Inertia, I, | 7.97 slug — J? | 

















C. DESCRIPTION OF BLUEBIRD 


The Bluebird aircraft was acquired as a test bed for guidance and navigation 
systems. Ultimately, these systems will be installed on the Archytas aircraft. The 
Bluebird is a conventional aircraft that will be used to test systems for a similar 
configuration to the Archytas in forward flight. The aircraft model was developed 
in the same manner as for AROD, as will be described in Chapter III. Physical 


characteristics for the Bluebird are given in Table 2.3. 


Ill. AIRCRAFT EQUATIONS OF MOTION 


A. NOTATION 


In this section the notation used in modeling the equations of motion is intro- 
duced. This notation is common in the field of robotics (see [Ref. Sil 91] and [Ref. 


Cra 86]). and is shown below in Figure 3.1. The following conventions are used: 





Figure 3.1: Relative Position of Coordinate Systems 


{A} represents the coordinate system with basis vectors, r4,ya, and z,. 
e “Po represents the position of point Q, expressed in {A}. 
e “Vo represents the velocity of point Q, measured in {A} and expressed in {A}. 


e °(4Vo) represents the velocity of point Q, measured in {A}, and expressed in 


{ B}. 


e 2 is a rotation matrix, also called a direction cosine matrix. A free vector 


in one coordinate system, that is a vector that can be moved parallel to itself 
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without change e.g., Vo can be expressed in another coordinate system by 


using the rotation matrix: 
“(CVo) = BR (PVo) 
e “Qp is the angular velocity of the {B} coordinate system with respect to {A}, 
and expressed in {A}. 


e °(40 8) is the angular velocity of {B}, with respect to {A}, and expressed in 


1B}. 


e Additional information can be added to the subscripts i.e., 4 Pgo is the position 


of the origin of {B}, expressed in {A}. 


B. COORDINATE SYSTEMS 


In order to derive equations of motion for a rigid airplane, the following coor- 


dinate systems and assumptions will be used: 


e {U} represents the inertial tangent plane coordinate system attached to Earth. 


{B} represents the body fixed coordinate system. 


e All sensors are located at the c.g. (This assumption will be lifted 
in a later section) 

e The mass of the aircraft remains constant. 

e Given a vector v, its derivative with respect to {B} is denoted as £(-) 
and 


its derivative with respect to {U'} is denoted as (') 
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The {B} coordinate system is defined with Xg as the thrust axis. A positive roll 


rate, p, is clockwise when looking in the positive X direction. The positive direction 


for Yg, the pitch axis, is out the right wing . A positive pitch rate, q, is defined as 


clockwise looking in the positive Y direction. The Zg axis is the yaw axis, and a 


positive yaw rate, r, is defined as clockwise, looking in the positive Z direction. 


To simplify the notation in places where it becomes cumbersome, The following 


definitions are introduced: 


c. 


@ vg represents the velocity of an arbitrary point, Q, measured and expressed in 


{U}. 


® vgo represents the velocity of the origin of {B}, measured and expressed in 


ee } 1.€. UVeo = UBO. 


e ’vug represents the velocity of point Q, measured in {U} and expressed in {B}, 


1.e. B(UV)) = avo: 


® wg represents the angular velocity of {B}, measured and expressed in {U}. i.e. 


UN, = wp. 


e 'we represents the angular velocity of {B} measured in {U}, and expressed in 


Wea.e. “(- Nz) = we. 


SPATIAL ORIENTATION 
1. Euler Angles 


The spatial orientation of a rigid body [Ref. Ju 92] can be defined by the 


three Euler angles, ®,9, and W called roll, pitch and yaw and defined in Figure 3.2. 


The Euler angles,in turn, can be used to define a rotation between two coordinate 


1] 





Figure 3.2: Z-Y—X Euler Angle Rotation Sequence 
systems. This rotation is obtained using Euler’s theorem: 


Any number of rotations about different axes through a point must, in 


the end, remain equivalent to a single rotation. 


For the case of conventional aircraft, a 3-2-1 rotation sequence is used [Ref. Sch 92], 
where the aircraft 1s yawed, pitched, then rolled. In the cases investigated here, © is 
small, and in steady state flight is equal to the angle of attack, a. @ can be expected 
to be anywhere from + 60deg in normal flight and can be anywhere from + 180deg in 
acrobatic flight. WV represents the heading of the aircraft and of course can range from 
0 to 360 deg. The transformation from inertial coordinates{U}, to body coordinates 


{B}, is carried out as follows, and is shown in Figure 3.2. 
1. The inertial coordinate system is represented by the vector “V, with the com- 
ponents x, y, and z. The first rotation is made about the z axis through an angle 
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Cs 


W. Now the vector is expressed as *V with the components r2, yz, and 22. Since 
the rotation was about the z axis, the z2 component remains unchanged. The 


resulting elemental matrix is: 


cosW sinWV 0 
M(W)= | -sinW cosW 0]. (3.1) 
0 eae! 


2. Now the rotation is made about the new y axis, y2, through an angle O. This 
results in a third coordinate system with the vector expressed as °V, and having 
components 23,43, and 23. This rotation does not change the y3 component. 


The resulting elemental matrix 1s: 


cosO 0 —sin9O 
M(O)= 0 ] 0 (3.2) 
sinQO 0 cosO 


3. Lastly, the rotation about the 73 axis through an angle ® is made to give the 
vector expressed in body coordinates, ?V. Now the resulting matrix is 
l 0 0 
M(®)=1|0 cos® sin® |. (3e3)} 
0 -—sin® cos® 
When the matrices are multiplied together in the correct sequence, M(®)M(O)M(¥), 
the result is the BR direction cosine matrix, expressed in terms of Euler angles as 


shown 


cosW cosO sin V cosO —sinO 
cosW sinOsin® — sinVcos@ sinOsin®sinW + cosWcos@ cosOsin® |. (3.4) 
cosW sinO cos® + sinW sin® sinOcos®sinV — cosWsin® cosOcos® 


The next step is to develop the kinematic differential equations that de- 
scribe the change in Euler angles with time. Following the method used in [Ref. 
Sch 92], the matrix of differential equations, 2, can be written as a sum of individual 


Euler angle rates: 


£o 0 0 
Q = M() “0 + M(%)M(O)| £0 | + M(®)M(O)M(W)| 0 |. (3.5) 
0 0 {y 
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When the matrix algebra in Equation 3.5 is done, the resulting kinematic differential 


equations for p,g,and r are given as: 


q}=1!10 cos® cosOsin® © (3.6) 
r 0 —sin® cos® cos® Ni 


The matrix on the right hand side of Equation 3.6 is invertible for all O 4 7/2, and 


can be used to solve for the Euler angle rates, ,0 and WU: 


o 1 sin@tanO cos® tanO Pp 
QO} =]| 0 cos & —sin ® q |. (3.7) 
y 0 sin®secQO cos® secO r 


By integrating Equation 3.7, the time history of the Euler angles can be obtained. 

The Euler angle method has one drawback. In the kinematic differential 
equations derived, a singularity occurs for some particular value of either ®,0, or W. 
In Equation 3.7, this singularity occurs at 0 = +7/2, which means the that the 
coordinate transformation is useful for an aircraft in horizontal flight, but useless 
for an aircraft which requires a vertical take-off or extended periods of vertical flight. 
Thus, another type of transformation is necessary for aircraft that spend considerable 
time flying at conditions near 0 = 7/2. 

The first alternative is simply to use another one of the 12 possible Euler 
angle transformations. It has been shown that the rotation matrix, R, is made up 
of sequential rotations and can be characterized as the product of three individual 
matrices where 


R= M,(03) Mg(@2) M.(41), (3.8) 
where the rotation sequence (a, 2, y) represents one of the combinations of integers 


(1,2, 3)e( ce eo ee ah 


(2, 1,3); (2,3. 1), Quo nee) 





and the A/,(0@;) are the rotation matrices 


l 0 0 
M,(@) = | 0 cos@ sin@ 
0 -—sin@ cos8é 
cos? 0 —sin@ 
M,(@) = 0 1 0 (3.9) 
sinf? 0 cos@ 
cos? sin@ 0 
M3(@) = __ io 0 
] 
One can see that by substituting in (3,2,1), and ®©,0, and W for 6), 02, and 63 re- 


spectively, the matrix shown in Equation 3.4 is obtained. Euler angle transformation 
matrices for each of these combinations have been calculated and are tabulated in 
[Ref. Ka 83]. The best Euler angle rotation sequence for an aircraft with flight at 
© = 7/2 was selected as a 2-3-1, or Y-Z—X, sequence. This sequence will allow flight 
at O = 7/2 with no corresponding singularity in the kinematic differential equations. 
Note that this matrix is invertible for W 4 7/2. The 2-3-1 rotation matnx, R, 1s 


described in terms of Euler angles as 


cos O cos V sin VU —sinOcos V 
—cosOsin Vcos®@+sin® sinO cosWcos® sinOsin V cos ® + sin ® cos O 
cos O sin V sin ® + cos ® sin @ —cosVsinO —sin@Osin Vsin® + cosOcos ® 





(3.10) 
The kinematic differential equations can be found in [Ref. Ka 83] as 
® 1 1 —cos@sinW sin@®@sin WV p 
at 0 cos ® —sin® q flu 
w on 0 sin®cosW cos®cos ¥ i 


These equations can now be integrated to find the time history of the Euler angles. 
2. Quaternions 
Another choice for the expression of a body’s spatial orientation is the 
use of quaternions. Quaternions eliminate the disadvantage of the singularity in the 


second rotation that is associated with the Euler angles. Quaternions have been in 


lo 


use for quite some time, having been discovered by Euler in a search for complex 


numbers [Ref. Mo 84]. A quaternion is defined as [Ref. Mo 84]: 
q = Bo + if + 82+ kBs, (Sa 


where the parameters are represented by various authors as S, a, b, c by [Ref. Ro 58], 
x, €&, 7, and ¢ by (Ref. Whi 59], and q4, 91. g2, and q3 by [Ref. Sil 91]. 

The components §o, 21, 82, and $3 are real numbers and the terms2,j, and k 
are defined in the typical manner for complex numbers, where 
ye =—)| = 
2 =) 
be ==) = 


The norm of a quaternion, g"q, is required to be 1: 
WV =¢ ¢=))45, +2 


since g* = By — 73; — 782 — k G3. 
It can be shown that A can be represented as follows: 


My. Py TF i13 

B 

git = || Toye pets (3.13) 
Pai” 1390" "ra3 


where 


= 2 22 pe 2 
les! re a ae eo 


re = 2(8)82 + 8o33) 
TASS. == 2( 8183 = Bo Bo) 
Oe A 1 B2 — Bo 83) 
[22> = ks a eh a SE _ Be : (3.14) 


I! 
to 


123 2( 3283 + BoP1) 


r31 = 2(8)83+ BoP) 
Ls ya 2( 6283 = Boh) 
rs33 = 83-6? — 83+ B3 


All that is required now is to determine the kinematic differential equations using the 


quaternions. 
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By expressing subsequent rotations from one coordinate system to another, 


where /' orients F; to F', and 8" orients F2 to F,, an algebraic approach [Ref. Mo 84] 


can be used to relate F, to F. 


R(B) 
R;;(B) 


R(B") R(B") 
Er (oitene | 


(3.15) 


Now the /’s can be expressed in terms of 3’’s and (”’s with the following result 


RB"), 

1D weriths eiie 
Bo B3 — By 
—-Ps Po fy 
ch aap es 


(3.16) 


(3.17) 


By regarding the second rotation in Equation 3.16as infinitesimal, the following result 


= 
where 
Bo 
By 
R — 
ye 
Bs 
is obtained 
B= 
where 
Bo 
_ | A = 
B a | p 
Bs 
and 
0 
le . 
W3 
Equation 3.18 can be rewritten as 
B= 


| R(ws)B 
Bo 
: Bnnal eo g — 
Bs 
Wy —W20 —W3 
0 Wa. 5 
meg Bie 0 Wy 
Wy —-Wy 0 
= R(B)ur, 


(3.18) 


| (3.19) 


be AES oO ce a 


(3.20) 


(3.21) 


and in this form can be integrated to give the time history of the orientation of the 


body. 


Ii 


Using quaternions has the following advantages over Euler angles in repre- 


senting spatial orientation of a rigid body: 
e Four states required to express the spatial orientation. 


e Requires almost 30 % fewer calculations [Ref. Ro 58], mainly because no non- 


linear, trigonometric equations need to be calculated. 


e No singularities in Equation 3.21 at any body attitude. 


D. DERIVATION OF EQUATIONS OF MOTION 

The derivation of equations of motion for a general six degree of freedom airplane 
model can be divided into two parts. The first part is simply the determination 
of the equations of motion for any rigid body in space. It is dependent only on 
the linear and angular momenta of the body. The second part is the calculation 
of aerodynamic, gravitational, and thrust forces on the airplane. These forces are 
particular to a certain aircraft and in general can be represented by the stability and 


contro! derivatives described later in the thesis. 


1. Linear Equations 
Equations for linear motion can be calculated using Newton’s law, F' = ma. 
Since most of the sensor information available for feedback to the control and nav- 
igation systems is available in the {B} coordinate system, the terms for linear ac- 
celerations, as well as forces and moments, will be expressed in the body coordinate 
system. First the position of the aircraft c.g. is determined as ¥ Pg9. Then Coriolis’ 
theorem is applied to obtain linear velocities for the aircraft. Coriolis’ theorem is 


then reapplied to derive the expression for linear accelerations. Then 
U 4Suy 
VBo = Po. (322) 
Both sides of Equation 3.22 are premultiplied by PR to get: 
BR’ Veo = GR” Pro 


Or 


Bupo = ®Pgo (3.23) 


Now consider Coriolis’ theorem 


A=TAtuxA (3.24) 


where A and fA use the notation for derivatives previously defined in Section A. The 
term w x A represents the difference between the relative velocity as measured from 
the rotating and non-rotating axes [Ref. Gre 88]. 


Equation 3.24 is applied to ?vgo in Equation 3.23 to get: 


Figo = <P 0p0 + Pus x Popo. (025) 
Newton’s law can now be written as 
UF = ma 
= tid Uso. (3.26) 


We 


where 'F is the total external force applied to the aircraft c.g. Equation 3.26 is 


premultiplied by #R to obtain the result: 
a oh = m ase "bBO 
= mM EF rea) (3.250 
when ?%go is substituted into Equation 3.27, the final result for ?F is 


d 
Br = m (T° eBo + Bip x BBO) 


=> m = UBO +m Pup x Bp. (3.28) 
2. Angular Equations 


The equations for angular accelerations are derived using Euler's law for 
preservation of angular momentum. These equations are also derived for the aircraft 


c.g. by applying Coriolis’ theorem to the equation for Euler's law: 
sea Naor (3.29) 


where “Lgo is the angular momentum of the aircraft and “ Ngo is the total ex- 
ternal moment applied to the aircraft c.g. Euler’s law can be rewritten in {B} by 


premultiplying Equation 3.29 by 2 R to get 
7 Lao = FR Noo. (3.30) 
Using Coriolis’ theorem in Equation 3.24, BIpo can be rewritten as 


d 
© ae re =, LBo =F aR: x ibe). (3250 


The angular momentum, ?Lgo, is defined as the product of the inertia tensor, Jp, 


and the body’s angular velocity, ?wg, or 
A 
BL = Ip Pwe + Ip Owe, 
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where Jp and "wr are the moment of inertia and the angular velocity of the propeller, 


respectively. When this term is substituted into Equation 3.31, the result is 


; d 
* Lo = = la we + Ip? wr) a sie x (Ip? wp + Ip? wp), (3:32) 


where /g is the inertia tensor for the aircraft and /p is the inertia tensor for any 
significant spinning object on the aircraft, such as a propeller, turbine , or other 


B 


rotating disk. The term, “wp, is the angular velocity of the rotor, expressed in {B}. 


We can carry out the differentiation in Equation 3.32 to get 


: d d 
“Lo = Ip "wp + Ina wr ST iy ap (Ip? wp + Ip? wp). (3°33) 


However, since d/dt(?wg) = 2g and d/dt(?wp) is assumed to be very small, Equa- 


tion 3.33 results in 

Bhgo = IpPwp + Bwp x (Ip? we + Ip?wr) (3.34) 
Now the result in Equation 3.34 can be substituted into Equation 3.29: 

® Ngo = Ip?wg t+ Pw X (Ip? wg t+ Ip? wr). (3.35) 


The term Jp?wpr can be disregarded if it is insignificant compared to Ig and PD 
[Ref. Ros 79]. 
3. State Equations 
In the preceeding sections, kinematic equations for the motion of a rigid 
body were derived in matrix form. These equations can be assembled into a state 
space representation of the kinematic equations of motion. First, Equations 3.28 and 


3.35 can be written as 


| oF | 7 | m +," upo +m ("wp X “vBo) (3.36) 


EIN) Tp? wp + Bop x (Ip? wp + IR? wR) 


Al 


Equation 3.36 can be rearranged to yield 


= m (—P wp * 2 Bo) +7 F) 


Le a (3.37) 
dt | Ip?wp ~ | —Bup x (Ip?we + IpBwr) +7N) | - 


The terms on the left hand side of Equation 3.37 can be normalized by multiplying 
by 1/m and Jj’. with the final result of 


d | PuBo | = | —Buyp x Bupo 


2 
dt| Fw lo se (ie Ip? i 
B —Ip “wp x (Ip? we + [Ro wr 


_ 
) ze f= BN | : (3.38) 


E. EXTERNAL FORCES AND MOMENTS 

Section D. gives the equations for kinematic motion, as shown in Equation 3.38, 
for any rigid body. Now it is necessary to distinguish between the different platforms 
to be modeled in order to give an accurate representation of the aircraft. This is 


~ 


achieved by computing the actual ?F and ?N acting on the aircraft. These forces 
and moments are those due to gravitational, propulsive. and aerodynamic effects, 


written as 


| we | = | ° Forav + °F prop + °F 4ERo (3.39) 


BN ® Nprop + ™ Nero 

1. Aerodynamic Forces and Moments 
The aerodynamic force and moment terms are determined by using a first 
order Taylor series expansion around a given nominal operating point. This operating 
point is the state chosen to represent the aircraft’s flight condition. Each term in the 
series is a partial derivative of ?F or ?.N with respect to the aerodynamic variables, 


u/U, a, B, p, q, r (Ref. Sch 92, Th 89]: 
WV on 6Fue +é6Par 0), hea (3.40) 
Similarly, moment terms can be written as 


NaERO = 6 Niwa’ + 6 Nix! + 6NxA+ No, (3.41) 


LS) 
tro 


where z’ is given by 


7 ila pe ge 6b ) 
t= 18% oF oe 20 ee, 
and zr’ is given as 
r’ = (3, d). (3.43) 
Control inputs are represented by the vector A: 
(N= Gono non| (3.44) 
where 6, 6,, and 6, are the elevator, rudder, and aileron inputs, respectively. 
Equations 3.40-3.44 can now be combined as follows: 
Aa OC OC 
| Ww, = 95(50 ste tat + 2,4 + Cro}, (3.45) 


where 9 = 1/2pV?, S = diag{S,S,S,Sb,Sc,Sb}, and C is the matrix of non- 
dimensional stability derivatives differentiated with respect to the terms defined in 


Equation 3.42, 3.43, or 3.44. 2 is defined as: 


CONC. iCamce Cy 
Cyre fey. «Cy. Cy. = Cy 
Crewe De ODO n,, “CH, 
Gracie cin Ser. Cr, Ci, 
OMCs Cas Cn, Cn. Cm, 
oF. om C5: One Cr, Cy 


ge is very similar to £5 ge except that only the a and B terms are normally computed, 


leaving a 6 x 2 matrix rather than a square matrix. The term 2 is defined as 


Cie CL © Ci, 
Gye Cy Cy, 


Cio Cin (Soy, 


Cis OF One 
CoC, 


on Cra Crs 
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C'ro 1s defined to be the vector of steady state coefficients: 


C'po 
Cyo 


_ | Cro 
Cro — Ga ’ 


C'n0 
Cro |! 


representing conditions in trimmed, balanced flight. This definition is similar to the 
definition used by Roskam [Ref. Ros 79]. In other references, the term Cro can 
refer to the nominal value of the coefficient at a = 0. However, in the Taylor series 
expansion it is more natural to use the first definition of Cro; therefore, it will be 
used in the following derivation and modeling. The stability and control derivatives 
are usually computed in the so-called wind axis coordinate system. The wind axis 
coordinate system, {W’}, is defined as the coordinate system that results when the 
Ip axis is aligned into the relative wind. This axis will not be aligned with the body 
coordinate system since the aircraft flies with an angle of attack, a, and can have a 
sideslip angle, 8. The transformation from {W} to {B} is performed in the same 
fashion as the Euler angle transformations mentioned earlier. The rotation matrix, 


WR, is a function of a and 8, and is expressed as 


cosacos3 ~cosasin? —sina 
R= sin 3 cos 3 0 (3.46) 
sinacos3? —~—sinasin3 cosa 


The rotation from {W} to {B} is made by premultiplying the force or moment vector 
by 3,R. Additionally, since the lift and drag are defined as positive along the negative 


zp and rpg axes, we define F'agro and Nagro as 


—D l 
F’ageRo = Ve and NAERO = ™m : (3.47) 
—L n 


In order to write Equation 3.45 in state space form, state variables must 


be defined. The most commonly used notation to use for the state vector 1s to use 


(3.48) 


| 


However, the terms z’ and z’ in Equations 3.40 and 3.41 cannot be used directly as 


states. Define 


M:2-72' 


M:22' (3.49) 


where 


M' > diag{1/Vr, 1/Vr, L/Vr, b/2Vr,c/2Vr, b/2Vr } 
and 


wa | Hen 0 Muto s 
~l0 0 d/(2Vr) 0 0 0 


are matrices of appropriate dimensions. The complete expression for ?Fagro and 


BN aero can now be written as 


BF aero _ -@ was 0 OC ! OC "rhe OC 
| anit | =a5| 0 BR (Cr ta iM t+ a Mat 2,4} (3.50) 


and can be substituted into Equation 3.38. 
2. Other Forces and Moments 
In addition to the forces and moments due the aerodynamics of the aircraft, 


forces and moments due to the propulsive and gravitational forces must be considered. 


Gravitational forces acting on the aircraft, Fray, can be found by premultiplying 
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’ Forav by the appropriate rotation matrix, where 


0 
UForav =| 0 
mg 
Then 
” Forav = DRY Forav- (3-309 


Propulsive forces and moments, ? Fprop and ? Npgop. are computed directly in {B} 


and can be expressed as: 


iy 
? Fprop = | Ty (3.52) 
a7 


and 
> Nprop = | Tm |, (3.58) 


where 7;’s represent the forces or moments due to thrust. Computation of propulsive 
forces and moments depends on each particular engine installation, and must be 
determined for the individual aircraft modeled. 


Equations 3.51, 2., and 2. can now be substituted into Equation 3.38: 


d UBO = — FU px 0 2 uBo a 
dt Bop - 0 —8 17) (Fup x (2 ieweoe es TRPwr)) Bop 


where 
ae ° Forav 4 ° Fprop a 
ae 0 ® Nprop 
a 0 Ve f fle aC id 
eh ip | 95 { Crot Soare + Sear'e + aa } hh (S/gen 


In order to write Equation 3.38 in state space form, the terms associated with zr’ must 


be collected and moved to the left hand side, along with the other time derivative 
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terms, Pigo and Pw. Let 


2 i m OQ 
B _ W 


then the complete non-linear equations of motion for any aircraft can expressed in 


state space form as follows [Ref. Th 89]: 


d Buso er ~ Fu ex 0 4 
di| Pup | * 0 —? I5*(Pwe x (PIp?we + IR? wr)) 
B B 
M718, TGS 35M’ | | Be | +Mz"{ | ia | + 
B 


0 
FpRoP 
E Li epereye | é7 + BT GS(Cr, + aL A) | (3.56) 


" Pao = OR go, aean)| 
and 
A = S(A)P wp. (3.58) 
where 
te BTSs Mt (3.59) 


P is the position vector of the aircraft, and A is the matrix of kinematic differential 


equations based on Euler angles or quaternions. 
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IV. COMPUTER MODELING O& TG 
AIRCRAFT EQUATIONS OF MOTION 


A. IMPLEMENTATION OF EQUATIONS 
1. Procedure 

The AROD/Archytas and the Bluebird aircraft models require different 
non-linear equations of motion. This difference is due to the unique nature of each 
aircraft. In the case of the AROD, the angular momentum of the propeller is a major 
factor to be considered, while the aerodynamics effects in a hover are negligible. 
The Bluebird on the other hand is a conventional aircraft, and exhibits the opposite 
characteristics. Angular momentum from the propeller is small and, the aircraft 
requires the stability and control derivatives associated with aerodynamic flight. All 
the equations were implemented in a systematic manner following the same general 
approach for either aircraft model. The commercial products MATLAB and SIMULINK, 
©) 1990-1992 the Mathworks, were chosen for the modeling’, mainly due to the ease 


of expressing matrix equations. The model was constructed using the following steps. 
e Kinematic equations of motion were coded. 


e Gravitational forces, with the direction cosine matrix represented by Euler an- 


gles were added to the model 


e Stability and control derivatives were included in the model, as well as engine 


thrust, as appropriate for the aircraft being modeled. 


1A]l code is listed in APPENDIX C. 


bo 
CP 


e Engine and actuators were added to the model. 
e Sensors for control systems were added to the model. 


In the first stage of modeling, analytic linearization was also carried out to 
verify the computer calculations. This was done by analytically linearizing the matrix 
formed from the six non-linear equations governing the kinematic motion of the body. 
Nominal values were substituted in, and the results compared to the trimmed and 
linearized values obtained from the SIMULINK program’. The next step required the 
addition of gravitational terms, and analytic linearization was still manageable. The 
linearized matrix now included nine equations and nine states with the Euler angle 
direction cosine matrix, DCM, and ten equations with ten states for the quaternion 
DCM. Nominal values for each case were again substituted into the linearized ma- 
trix and compared to the linearized model derived from SIMULINK. The inclusion of 
the stability and control derivatives and the thrust terms presented a problem that 
was much too cumbersome to linearize analytically. Verification of the data at this 
stage was accomplished by direct comparison of the dimensional derivatives result- 
ing from numerical linearization of the plant. In the case of the AROD the results 
were compared with data published by Sandia Labs [Ref. Wh 91]; for the Bluebird, 
eigenvalues were computed and then compared to eigenvalues obtained by classical 
analytic methods [Ref. Sch 92] and [Ref. Ros 79]. 

2. AROD Equations 

The AROD differs from the Bluebird primarily in that the aerodynamic 
forces and moments are negligible while the craft is hovering. The powered lift does 
present some special problems, the predominant difficulty being the gyroscopic cou- 


pling of the spinning propeller. Another important consideration is the moment due 
Data is tabulated in APPENDIX B. 


BS 


to the swirl effect of the air from the propeller striking the vanes mounted aft of the 
propeller. These forces and moments are computed, then substituted in for ?F and 
PN in Equation 3.38. 
a. Control Forces and Moments 

All of the applied forces and moments in the AROD are due to the 
four control vanes mounted aft of the propeller. The vanes can be moved in different 
combinations, as discussed in later in this chapter, in order to maneuver the AROD. 
The forces and moments, ?F and ?N, acting on the AROD are computed from a 


Taylor series expansion around a nominal hover point. Since aerodynamic terms are 


negligible, 
[anc |= aangga (42) 
where 
C= T0ew On onl 
eg; = 1/2pV, 
e V; is the induced velocity through the propeller [Ref. Pro 90] and 


Ve = T/(2Ap). 


A is the inlet area, where A = 3.14 f?. 
e R is the propeller radius, where R = 1.0 f°. 


Notice no aerodynamic forces act on AROD due to the movement of the elevator, rud- 
der, or aileron controls. Again, this is because the model of interest is only designed 
to fly in a stable hover. The other forces and moments involved in these calculations 


are due to gravitational and propulsive action, and are described next. 
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b. Additional Forces and Moments 
The gravitational force term, ? Foray, is computed in the {B} coor- 
dinate system with the rotation matrix, R. For the AROD, §R is determined by 


the 2-3-1 Euler angle rotation sequence defined in Equation 3.10 and is written as 


0 
?Forav =GR’Forav where”’Forav =| 0 |, (4.2) 
mg 
where 
—sinO cos WV 
BForav = mg} sin® sinV cos@+4+sin® cosO |}. (4.3) 


—sinO sinWV sin®+cosO cos ® 


The forces and moments due to the propeller thrust were determined 
experimentally [Ref. Sto 93] and are discussed in a later section. For the AROD, the 


propulsive force, B Fprop, is acting completely along the zg axis 


1S 
2 F prop = 0 : (4.4) 
0 


where F’r, is the total thrust, determined as a function of rpm. 
The moment resulting from thrust is ?Nprop and is given by the 


vector equation 
lr 
5 N prop =~) 0 (4.5) 
0 


where /7 is the rolling moment due to thrust. This rolling moment is due to the 
swirl of the air as it leaves the propeller and strikes the control vanes. This term was 
determined experimentally as a function of thrust, and is also discussed further in a 
subsequent section. 

When all the terms are collected, the total force and moment applied 
to AROD can be expressed as 


ie ae male 
aa Oe PROP GRAV 
fan ]={ 74886+| expen |t[ “o" | fe 


3] 


The force and moment terms in Equation 4.6 can be substituted into 


Equation 3.38 with the resulting state space equation given by 
d = UBO _ — ie x Bozo i 
FT a): —5- Fwe x (Ip? we + In?wp) 


I/m 0 ee ° Fprop ? FGRav 
| 0 jee | ( gGARZ.& + | Nee + 0 : (4.7) 


Equation 4.7 can now be written in the state space form and programmed using 


SIMULINK 


ad Bupo = Baye 0 *uBo 
dt Bup 7 0 —* 15 (Fwp x (F Ip wp > Ip? wr) Bop 


+ 
(ones [2am [ "2 J 


° NpRop 
3. Model Validation 
a. Kinematic Equations 
As discussed in the beginning of the chapter, the first step in the 


AROD model validation was to linearize the non-linear equations 


dg 


yj UBO = 1/m(—"?wp x ?vgo) (4.9) 
d 
= ws = ? Ip (—"we x (FIR? wr + “ Tpe ony (4.10) 


with the result given by?: 


d dv | | —wox Vo X bv 
dt | dw | 0 —FI5'(FIRewr x —(wox)? Ip t+ PIpwox) | | dw 


+ | : | dbwr. (4.11) 


—Wax ai ire 


The nominal values for the AROD hover operating point, 


10 
0 
0 
0 
0 


[0 | 


3See APPENDIX A for a description of the linearization process 
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can be substituted into Equation 4.11 to get 


000 0 0 0 
000 0 0 10 

d[év)_|000-10 0 0 bv 

al a | 000 0 0 0 _ Si 
000 0 O -1.615 
1000 0 1606 0 


where the linearized matrix associated with dwp is zero, since —wp x Ip = 0. These 
values match the values obtained by linearizing the model using SIMULINK. 

The inertial cross coupling from the propeller is evident from the lin- 
earized results. The values of —1.6152 sec”! and 1.6062 sec”! are defined as pitching 
moment due to yaw rate, m,, and yawing moment due to pitch rate, n,, respectively. 
The gyroscopic coupling is demonstrated by putting a step moment input of 1 second 
duration along the y axis and observing the time history of the angular rates q and 
r as shown in Figure 4.1. The AROD has a tendency to spin in the axis orthogonal 


to the input torque as shown by the motion r about the z axis. 
0.8 
0.6 
0.4 
02 


0 


p tadians/second and Phi, radians 





Time, seconds 


Figure 4.1: Gyroscopic Motion of AROD 
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b. Gravitational Forces 
The next step in modeling the AROD was to include gravitational 
effects into the model. This required the expression for ?Foray, given by Equa- 
tion 4.3, to be coded into a MATLAB function. This function block was then added 


to the SIMULINK model, as shown in Figure 4.2. The equation to be analyzed at this 





Figure 4.2: Modeling of Gravitational Effects for AROD 


stage was 
B 
si Bupo — wp x Bape ae a 
Fh "wp | = | Big (—Bwex (FlaFwr t+ Pig? we) |; (4.14) 
A S(A)?wep 


where S(A) represents the kinematic differential equations, defined in Equation 3.11. 
Notice that since gravity acts through the c.g., no ?>Ng¢ray term exists. Next, Equa- 
tion 4.14 was linearized both analytically and with SIMULINK. The analytic result, 


given in state space form, was 


1 dv —Wy X Vo X f(A) dv 
°F COs = 0 PTE (BI pPwr x —(wox)?lg + ? IpwoX) 0 dw 
dA 0 G(A) h(A.w) |} 6A 
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0 
ae to eae | dwp. (4.15) 


0 





In Equation 4.15, the functions f(A),G(A), and h(A,w) are partial derivatives of 
Equations 4.3 and 3.11 with respect to A, or 0/OA. The function f(A) is the linearized 


expression for ? Foray given by Equation 4.3, where 








9 BR 9 — sin Ocos V 
ley) = aA oe =IFq sin © sin V cos ® + sin ® cos O (4.16) 
a —sin Osin V sin ® + cosOcos ® 
or 
0 — cosO cosV sinQ sinW 


g| —sinOsinWV sin®+cos®cosO cosOsinWcos® —sin®@sinO  sinOcosV¥ cos® 
—sin® cosO — sinOsinY cos® —cos®sinO — cosOsinV sin® —sinOcosY sin® 
(4.17) 





where f(A) is evaluated at the nominal condition, Ag. The function 
G(A) = 0/0w(S(A)?wp 


is derived by linearizing the kinematic differential equations in Equation 3.11. The 
function h(A,w) = 0/0A(S(A)?wp) is not presented since for wo = 0, as in steady 


state cruise, h(A,w) = 0. The matrix G(A) is derived from the following equation 


9 1 —cos@tanWU  sin®@tanV p 
eA) = 5 (S(A) wa) =. 0 cos®secW —sin ®sec V q (4.18) 
0 sin ® cos ® r 


0 cos®secV —sin®sec V 
0 sin ® cos ® 


1 —cos®@tan V sin ® tan UV 
(4.19) 
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These matrices were assembled as in Equation 4.15 and were evaluated for the nominal 


flight condition. The values for zo for this flight condition were given by 


10 


i) 


lo = 


(4.20) 


So 2 SS 


fete 1009 || 


After substitution into Equation 4.15, the result for the linearized equations of motion 


was 

000 0 0 0 Oo oO 0 
000 0 0 10 0 0  —=esoeed 
(0 Oe=l0m 0 0. 105 325 

$v 000 0 0 i 0 0 0 

a =|000 0 O -1615 0 0 0 (4.21) 

5A 000 0 1606 0 060 0 0 
ie0i 0 le 0 oO 0 0 
0.50/20 50 oeteet 0 oO 0 0 
000 0 0 1 0 0 0 


This result was essentially indentical to the results obtained by using the trim and 
linearize functions in SIMULINK. The gravitational effects acting on the aircraft were 
examined by running a simulation of the non-linear model and comparing the results 
to those from Equation 4.11. The expectation is that a body falling in earth’s gravity 
will experience an acceleration of 32.174 f/s* and as shown in Figure 4.3, this expec- 
tation is realized in the model. In Figure 4.3, at the end of the 10 second simulation, 
the vertical velocity, u in this case, is ~ 320f/s, as was predicted. 
c. Additional Forces and Moments 
The last step in modeling the AROD was to include the forces and 


moments due to propulsive and contro] action that act on the aircraft. The data used 
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u. and w, teet/second 





Time, seconds 


Figure 4.3: Velocity of AROD, Gravitational Effects Included 


in modeling these control forces were collected experimentally [Ref. Sto 93] and then 
curve fitted to a function. The measurements of rolling moment, thrust, rpm, and 
vane position were taken for various configurations. The data collected were then 
reduced and the required characteristics were computed. The accuracy of the AROD 
model depended on very accurate modeling of thrust as a function of rpm, moment as 
a function of thrust produced, and moment produced by deflecting the control vanes 
in the different combinations. 

Thrust and rolling moments were measured directly at different power 
settings ranging from 3000 rpm to 7600 rpm, with a power setting of ~ 6400 rpm 
giving a thrust of ~ 90/b;. This thrust is aprroximately the force required to maintain 
a hover for the basic AROD configuration. Figure 4.4 shows the linear curve fit 
through the thrust vs. rpm data. The curve was fit using data from 5000 rpm to 
7600 rpm, as this was expected to approximate the normal operation range in flight. 
The thrust and moment data are plotted in Figure 4.5 along with the line fit through 


those data points. The best fit, by least squares, for the data in Figure 4.4 was given 
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T=0 0297*rpm - 1047 
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Figure 4.4: Thrust vs. RPM for AROD 
-56 
58 eos Moment due to Thrust -@-- 
; Least Squares Ft —— 
& 
6.2 
2 64 
= 
= 6.6 
ef 
© 
= 


6.8 


72 M=-0.0$42°T 0.914 


-7.4 





-7.6 
85 90 95 100 105 110 115 120 
Thrust, b 


Figure 4.5: Thrust and Moment Data for AROD 


by 


zt 


Fy, = 0.0297 6m — 104.7, (4.22) 


where 4,5. represents the rpm at a given throttle setting. The best fit for the data 


in Figure 4.5 was given by 
Nr = —0.0542 Fr, — 0.9138. (4.23) 


These equations were then used in modeling the engine response to a rpm setting. 

The engine itself was modeled using a second order transfer function 
from the servo position to 6;pm, based on the Sandia Labs [Ref. Wh 87] model. The 
transfer function for the engine was given as 


A 2 
lig oe 


+ — 4.24 
s? + 2Cw,ns + w? ‘i ee) 


rom 


where Kg = 900, w, = 5rad/sec, ¢ = 1.0, and é7 is the throttle servo position. Since 
the actual throttle position is set via a radio link and tests have not been set up to 
model the response, it was ignored in the model. 

The engine servo could be modeled and a transfer function from com- 
manded input to servo position was determined. This transfer function also was 
determined by Sandia Labs in the original AROD work [Ref. Wh 87] to be 


uw)? 


6; = > 4.25 
peg? 4 26 WS + w2 a ae) 


where w, = 20rad/sec, ¢ = 0.6, and ur represents the commanded input to the servo. 
It should be noted that Equation 4.25 was also used for the control vane servos, since 
all the servos in the AROD were identical. 

In order to model the moments due to control surface deflection, 
BN Control computation of control vane effectiveness was necessary. Again, the data 
gathered by [Ref. Sto 93] was used to compute vane effectiveness for the AROD 
model. Figure 4.6 shows the moment data for 75% thrust. This power setting cor- 
responds closely to the thrust required to maintain hovering flight. A line was fitted 


through the data points, and the slope of this line was used to determine the rolling 
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I=0.5523°V_AVG - 5.7248 
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Figure 4.6: Moment Due to Vane Deflection 


moment,/, for a given vane deflection. The equation for the curve was 
l= 0.5523 Vave — 5.7248 (4.26) 


for 6580 rpm, where Vayc was the average vane deflection, 1/4 5°?_, V;, in degrees. 
This averaging was required since the individual vanes were at slightly different posi- 
tions with respect to the commanded position. Computing the vane effectiveness for 


roll was done by differentiating with respect to the vane deflection, 


a sft — ly 
* OVave | sideg«S 





(4.27) 


The rolling moment should be non—dimensionalized for use in the equations of motion 


given by Equation 4.8. This was done by applying the following equation 


ls. P 
Cis, = 1/2pV2Sb (4.28) 


To use the measurable quantities available, redefining the terms in Equation 4.25 was 


necessary. The representative area, S, is defined as the inlet area to the duct, A. The 
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characteristic length, b, is defined as the propeller radius, R. The velocity term was 
defined as induced velocity, V;. With these terms the non-dimensional coefficient Ci,_ 


could be calculated as 


ls 
‘ie 4.29 
fa 1/2pV2AR ) 


and was computed to be C);, = 1.438/rad. The rolling moment was measured during 
the testing and was easily non—dimensionalized in a form suited to computer modeling. 
However, no measured data existed for the pitching and yawing moments. This 
required the estimation of pitch and yaw vane effectiveness by a simple ratio technique. 
First, it was assumed that the vane effectiveness, ls, for two vanes (2V), would be 


exactly twice the Js, for one vane, and 1/2 the Is, for four vanes, (4V). This would 








make 
al | eo 
— Yo | ie 4.30 
Vive 2 OVaveG Ja ses 
or numerically, 
Ol a 
, = 0.2762————— 4.31 
(Wave Jov = ra ( ) 
Now, the ratio of the dimensional derivatives to the moment arms was set up as 
al am 
(arava )2¥ _ (a¥yv) (4.32) 


l, ly 
where 1, was the distance from the c.g. along the x-axis, to the midspan of the control 
vane, a distance of 9.0 in. The distance J, was the distance from the c.g. along the 
y-axis, to the 1/4 chord of the control vane, a distance of 15.43 in. The calculation 


was performed and the result for the pitching moment was 





om ne = lbs 
7 (4.33) 


ah 47 35 
OVavG 


This was non—dimensionalized using Equation 4.28 with the result of C,,, = 1.2332rad7?. 


Moreover because of the symmetrical design of AROD, Cn, = C,,.. The results for 


4] 


vane effectiveness are important to a high fidelity model and are presented along 


with other relevant data in Table 4.1. The non-dimensional derivatives in Table 4.1 


TABLE 4.1: NON-DIMENSIONAL DERIVATIVE DATA FOR AROD 


[rad | deg | Positive Vane Deflection | 
LG, [1438 [0.01 | +++ al 


[Gna [1233] -002I5[ CSV | 
[Cn, [1233 [0015 SCO | 





were substituted into the term (OCr/0A) in Equation 4.7. Written as a matrix, the 
non-dimensional derivatives are 
) 
aCr Coy 0 


N= | fC, BU 
ee (0 Cs, 


m 


NH 
4 


(4.34) 


SN 
re) 


Equations 4.22, 4.23, and 4.34 were added as a separate block in the model, resulting 


in a block diagram shown in Figure 4.7. 





Figure 4.7: Non-Linear Equations of Motion Model 


In the SIMULINK linearization results, no change was expected in the 
A matrix since no aerodynamic terms were added to the model. Rolling motion is 
shown in Figure 4.8 as an example of the unstable nature of the AROD without a 


stability augmentation system. 
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Figure 4.8: Rolling Motion For Complete Non-Linear Equations 


4. Bluebird Equations of Motion 
a. Kinematic Equations 
Modeling the kinematic equations of motion for Bluebird was accom- 
plished using the same procedure as was used in the AROD modeling. A slight dif- 
ference arose in the derivation of the angular acceleration equation for Bluebird since 
the term, [pwr, for the angular momentum due to propeller rotation is negligible. 


Thus, the following equations were used in the first stage of modeling. 


d 
7 BO = —Puwp x Pugo (4.35) 
d 
= WB = 2ie\—"wp x PI p?wp). (4.36) 
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These equations were coded into a MATLAB function and placed into a block diagram, 


shown in Figure 4.9. Equations 4.35 and 4.36 were linearized analytically’ with the 





denv.m 


State _ dot 
lambda_dot 


Figure 4.9: Block Diagram of Kinematic Equations of Motion 
resulting state space equation 

da bv > = Win Uo X dv (4 37) 

dt | dw | 0 —B IF {—(wox)? Ig + P1g(wox)} dw |. 
This equation was evaluated at the nominal flight conditions, determined by the 


typical cruise condition for the Bluebird aircraft. A state vector for the nominal 


flight conditions is given by 


Ilo = 


(4.38) 


k____. 


4See APPENDIX A for the details 


fd 


These values of rg were substituted into Equation 4.37: 


0 0 0 0 —6.675 0 
0 0 O 6.675 0 —72.995 
d | dv 000 O 72.995 0 dv 
ol a | tO’  @ 0 0 a oy) 
0 0 O 0 0 0 
0 0 0 0 0 0 


These results were in complete agreement with the data obtained by trimming and 
linearizing the non-linear model with the SIMULINK program. The comparison be- 
tween Equation 4.39 and Equation 4.13 shows the absence of any angular rate cross 
coupling in Bluebird. The absence of the cross coupling terms results from the choice 
to ignore the angular momentum contribution from the propeller. 
b. Gravitational Forces 

After the basic kinematic equations Equation 4.35 and 4.36 were put 
into block diagram form, it was an easy matter to include additional blocks. The 
next block was a function block including the ? Foray terms. The model with the 
gravitational terms included is shown in Figure 4.10. The equations of motion to be 


modeled at this stage were 


d 1 
7 UBO = see ips x ag = eA (4.40) 
t m 
d 
7 WB = a rere x BIB? wp) (4.41) 
d 
Fi = S(A)Pwp (4.42) 


where S(A) is the set of kinematic differential equations based on the 3-2-1 Euler 
angle rotation in Equation 3.7. These equations were then linearized analytically, 


with the result 


d ov ae Vo X g f(A) bv 
az} be [=| 0 -Flg'{—(wox)? Ie + ?ip(wox)} 0 | | bw | (4.43) 
= 0 G(A) aon | LAN 
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Figure 4.10: Gravitational Forces Model 


where the terms f(A), G(A), and h(A,w) are based on the 3-2-1 Euler rotation se- 
quence. The linearization is done in the same manner as was shown in Equations 4.16, 


4.17, and 4.19 resulting in 


0 —cosO 0 
ie cosOcos@ —sinOsin® 0}, (4.44) 
—cosOcos® —sinOcos® OQ 
and 
1 sin@tanO cos®tanO 
Gv = 10 cos ® —sin® |. (4.45) 
0 sin®secQO cos®secO 
Note that 
O 
h(A, w) = 5A SUA) wale = 0, 
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since wy = 0. When Equation 4.43 is evaluated at the nominal condition zp given by 


72.9954 f/s 
0 
6.6757 f/s 
0 
fo = 0 ) 
0 
0 
0.0912 rad 
0 
the resulting equation is 
000 0 £4-6.675 0 0 —32.043 0 
0 0 0 6.675 0 —72.995 32.043 0 0 
000 O 72.995 0 0 —2.930 0 
d dv 000 0 0 0 0 0 0 dv 
—|édw}=]0 00 20 0 0 0 0 0 bw 
Ba) 5A 000 0 0 0 0 pw aN 
HOO 0 1 0 0.0915 0 0 0 
000 O 1 0 0 0 0 
000 0O 0 1.004 0 0 0 
(4.46) 


In this instance, as well, the results of the analytic linearization in Equation 4.46 
agree very closely with the computed results. 

A non-linear simulation of the system in Figure 4.10 should show an 
increasing downward velocity due to the gravitational effects of B Foray. One would 
also expect to see decreasing forward velocity due to the “drag—like” term that arises 
with the introduction of the angle, QO. This plot 1s shown in Figure 4.11. 

c. Aerodynamic Forces and Moments 

Completion of the Bluebird equations of motion model required the 
modeling of the aerodynamic forces and moments acting on the aircraft. A simple 
engine model was also developed. No analytic linearization was performed at this 
stage due to the increased complexity of the model. Verification of the computer 


results was accomplished by comparing the modes and eigenvalues of the computed 
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Figure 4.11: Gravitational Effects on Velocity 


plant with those resulting from substituting the stability and control derivatives into 
equations developed in [Ref. Sch 92]. 

The aerodynamic forces and moments as described in Equation 3.50 
were coded as a MATLAB function, then included as a block in the model shown in 
Figure 4.12. 

Next, it was necessary to premultiply all the blocks by y7', since 
this added the effects of the a derivatives. Now the important task was to calculate 
the stability and control derivatives using the general reference for the estimation of 
non-dimensional derivatives, DATCOM [Ref. USAF 60]. The stability and control 
derivatives were computed based the aircraft geometry and the control surfaces. These 
values are tabulated in Table 4.2 and in Table 4.3 where the non-dimensional force 1s 
listed on the left side, and the particular derivative is determined using the top row. 
For example, Cp, = 0.188 using Table 4.2. The terms in Cp, were also estimated to 
be C,, = 0.385 and Cp, = 0.03, with all other terms equal to zero since the aircraft 


is in straight and level flight at this trim condition. 
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Figure 4.12: Full Non-Linear Equations of Motion Model 


TABLE 4.2: NON-DIMENSIONAL STABILITY DERIVATIVES 





[cr [ol oosr[ 0 [0363 [0 | 0100 | 0 
[Cn [of 0 [ait 0 [am] 0 [470] 
[e. [ofoos7 | 0 [-oomr[ 0 [oom] 0 | 






The ? Fppop was based on estimating engine thrust for a 4 HP engine 


and a propeller efficency, np, of 0.65. Thrust, Jo was estimated using the equation 


_ 990 np HP pr 


V5 
Uo Po 


where py, is the density at the operating altitude, po is the sea level density, and Us is 
the velocity of the aircraft in f/s. The result was an engine thrust of To = 19.51b,, for 
a density ratio of 1. This could be directly factored into the equations of motion as 


Toor, where 67 was arbitrarily scaled from zero to one. 67 = 0 represents zero thrust 
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(4.47) 


TABLE 4.3: NON-DIMENSIONAL CONTROL DERIVATIVES 


a 
[Co| o0e5 [0 | 0 | 
[cy] 0 [oom] 0 | 
[c. [uml o [0] 
Pc [0 [0.0028 | 0265 | 
[c.[140f 0 | 0 | 
rc. | 0 [0.0329 [0.0377 | 










and é7 = 1 represents maximum thrust. 

The preceeding values were substituted into the appropriate MATLAB 
functions and the entire equations of motion model was trimmed, linearized, and 
then compared with the analytic results based on classical techniques for determining 
eigenvalues and eigenvectors. The eigenvalues are compared in Table 4.4. The results 
from the computed eigenvalues is very close to the eigenvalues derived by analytic 


methods. 


TABLE 4.4: COMPARISON OF EIGENVALUES FOR BLUEBIRD 


MODE COMPUTED ANALYTIC 





[LONGITUDINAL] 

[__ Phugoid | —0.010) + 0.4008) | 00473 20.4010) 
[Short Period | —3.9833 + 3.5521) | —4.008 = 3.5462) 
[_TATERAL [| 
[Dutch Roll | —0.538 4 3.0310) | 0.522 = 3.0191) 
[Short Period [3.5029 | 3.00 
[Spiral [0000.0 
















B. VALIDATION OF AN INDEPENDENT CASE 
Although verification of the model was accomplished at each stage, comparison 


of the results from the numerical linearization with linearized results from an inde- 
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pendent source was still necessary before the model could be considered completely 
reliable. The test case selected was the Cessna 172 documented in [Ref. Ros 79] as 
airplane A. The tabulated non-dimensional derivatives and given flight conditions 
were used as inputs to the non-linear model. The model was then trimmed at the 
specified flight condition of Vr = 219 f/s and Oo = 0. Using the resulting state and 
input vector, the model was linearized around thesed nominal conditions. With the 
linearized plant, eigenvalues were determined and compared to those tabulated for 
airplane A (see Table 4.5). Very little difference between the eigenvalues is seen in 
the longitudinal modes. Slightly greater differences are noticed when comparing the 


lateral modes, but these differences are still fairly small. Another very good method 


TABLE 4.5: EIGENVALUE COMPARISON FOR CESSNA 172 TEST 
CASE 












[Mode [Numerical [Independent] 
[Longitudmal [OT ~*™ 
[Short Period | —4130243005) | —413044300) | 
[___Phugoid | -0.0209 0.17947 | 0.02002 £ 0.1797) | 
[Lateral-Directional| | Cid 
[Dutch Roll | —0.0047 3.3080) | —0.S5NE3.300)_| 
[Spiral if daa09 
[Roll Response | 0.0109 (0.01005 _—| 


for comparing the results of the numerical linearization with Roskam’s tabulated data 
is to form plant and control matrices for the test aircraft, using the linear algebraic 
method taught in AE 3340 [Ref. Sch 92]. The resulting A and B matrices were then 
compared by applying step elevator, rudder, and aileron inputs and plotting the re- 
sults. The results for the step elevator input are given in Figure 4.13. These plots 
show very little difference in the vertical velocity, w. Differing amplitudes are shown 


for the horizontal velocity, but since the non-dimensional velocity, u/Vr, from the 
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dimensional derivatives was scaled to be equivalent to the state, u, computed in the 
numerical linearization, these magnitude errors are not indicative of a poor model, 
only a slight difference in the computed damping is shown. The results from the 
analytic model were scaled up by the nominal airspeed, Vr to compare with the nu- 
merically linearized results. This would have the effect of magnifying any differences 
between the analytic model and the numerical model. The natural frequency of both 


the analytic and numerical results are quite similar. 
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Figure 4.13: Comparison of Longitudinal Responses to Step Elevator Input 


Lateral responses to step rudder and step aileron inputs are shown in Fig- 
ures 4.14 and 4.15. Very little difference between the models is visible in these plots. 
The errors are shown in Figures 4.16 - 4.18 

It can be concluded that the results from the numerical linearization are quite 
close and furthermore that the equations used in the model are indeed correct. More- 
over, the linearized results presented for both the AROD and Bluebird aircraft should 
be considered accurate and are suitable for the Guidance, Control, and Navigation 


system designs that will follow. 
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Figure 4.14: Comparison 
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Figure 4.15: Comparison of Lateral Responses to Step Aileron Input 
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Figure 4.16: Difference in Analytic and Numerical Results, Step Elevator 
Input 
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Figure 4.17: Difference in Analytic and Numerical Results, Step Rudder 
Input 
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Figure 4.18: Difference in Analytic and Numerical Results, Step Aileron 
Input 
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V. SENSOR AND ACTUATOR MODELING 


The AROD engine and actuators were modeled as a second order transfer func- 
tions, based on data collected by Sandia Labs {[Ref. Wh 87]. Sensors were also mod- 
eled as second order transfer function based on data supplied by Watson Industries 
(Ref. WAT 93]. 

The complete non-linear model for an aircraft should include models of the 
sensors on board for measurement of acceleration, angular rates, pitch and roll an- 
gles, and headings. The inertial device chosen for the AROD project was the Watson 
Industries IMU-600D inertial measuring unit. This device contains a triaxial ac- 
celerometer, a triaxial rate sensor, two liquid pendulous devices (for bank and pitch 
angle), and a magnetic heading indicator. The characteristics of these devices must 
be accurately modeled, since it is the sensor output that drives the control system. 
In the rest of this section, the accelerometers, rate gyros, and inclinometers will be 


modeled, as well as the sources of error inherent in the design of the sensor devices. 


A. ACCELEROMETER MODELING 
The term accelerometer is not entirely accurate, since the device does not mea- 
sure true acceleration, but rather the difference between acceleration and gravity 
[Ref. Bro 64]. This effect is referred to as the Einstein Uncertainty Principle, and is 
represented in equation form as 
f=g-a (5.1) 
where g and a are the specific forces of gravity and acceleration of the aircraft and 
are measured positive downwards. The accelerometer model relevant to this equation 


is shown in Figure 5.1. The tri-axial accelerometer of the IMU-600D can be modeled 
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Figure 5.1: Typical Accelerometer Model 


as three simple single-axis accelerometers, as has been established through conversa- 
tions with the manufacturer [Ref. WAT 93]. A schematic representation basic device 


pictured in Figure 5.1 is modeled by an ordinary differential equation [Ref. Sil 91] 


ite 2 fi + =, = —v (5.2) 
where x denotes the displacement of the mass from its equilibrium position and 
y = g — ais the projection along the case axis of the vector sum of gravity and 
acceleration. The terms, 8,4, and m represent the damping, spring coefficient, and 
mass, respectively, of the device. The accelerometer described in Equation 5.2 can 
be modeled as a second order low pass filter, but the actual accelerometer has a flat 
response up to 1000Hz, so it was not modeled. A third order Chebychev anti-aliasing 
filter with a cut-off frequency of 20 Hz was added to the accelerometers. This filter 
is the device that was modeled. The Chebychev filter gives the advantage of a flat 


passband, and a very sharp drop off at the cut-off frequency. The equation used to 


o¢ 


describe a Chebychev filter [Ref. St 88] is 


i 


i? 2 ee 
\Hip(jw)| 1+ C2 (w)’ 


(5.3) 


where Cy is the Ni, order Chebychev polynomial, € is the parameter thats sets the 
ripple in the passband, and |Hzp(jw)|? is the magnitude of the filter. It was not 
necessary to compute the filter analytically, as SIMULINK provides a block function 
which performs the required steps based on the passband ripple of 0.ldb and the 
desired cutoff frequency of 20 Hz. The block diagram for the accelerometer model 
is shown in Figure 5.2. Included in this diagram are the blocks containing the error 
modeling, as well as the blocks used to correct for the Einstein uncertainty. Figure 5.3 
shows the linear synthesis model used to generate the accelerations to drive the sensor 


models. The synthesis model was derived from the Bluebird model discussed in 


Chapter IV. 





Figure 5.2: Accelerometer Modeling 
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Figure 5.3: Synthesis Model for Accelerometers 


1. Error Model 
No matter how well the sensor device is constructed, any accelerometer 
is subject to certain errors in the linear acceleration measurements. These errors 
can occur for several reasons; some as mechanical and others are due to the physical 
placement of the accelerometer on the aircraft. The mechanical errors accounted for 


in the IMU-600D tri—axial accelerometers are 


e Acceleration Bias 


— average readings for +lg and -lg loads 


e Acceleration Scale Factor error 


—average difference between readings from +lg and -lg loads 
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e Cross Axis Sensitivity errors 
—measurements due to the misalignment of an accelerometer with the appropri- 
ate axis 


—~measurements of off-axis accelerations are measured 


e Acceleration Noise Floor 


-threshold below which measurements of acceleration can not be made 


Errors in the measured accelerations can also occur if the accelerometers 
are not located at the c.g. of the aircraft, since arbitrarily located points on a body 
will experience additional accelerations due the the angular momentum of the body. 
A mathematical correction for the error will be presented later, and will describe the 
angular and linear accelerations of an arbitrary point on a rigid body. However, the 
correction will not be applied to the models here, since the correction is expected to 
have a very small effect on the sensor measurements due to the small displacements 
away from the aircraft c.g. 

Error terms are quantified in terms of full-scale measurements. The me- 
chanical errors are tabulated in Table 5.1 along with other important characteristics. 


Figure 5.4 shows the block diagram with error inputs applied to the accelerations 


TABLE 5.1: ACCELEROMETER CHARACTERISTICS 


[Acceleration Range [| -2ags—*d 
[ Acceleration Bandwidth | ___20 Hz 
[Acceleration Bias | 0.2% of Full Scale | 
[ Acceleration Scale Factor | 0.2% of Fall Scale | 
[ Acceleration Noise Floor | 0.0005 q's __| 
[_ Gross Axis Sensitivity | 0.5% of Full Scale | 
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Figure 5.4: Error Model for Accelerometers 


measured by the accelerometers. The output from the error computations in Fig- 
ure 5.4 is the measured acceleration from the accelerometer output to the control 


system. The cross axis terms are determined through a matrix 


; OL ees 
Tae eon) & (5.4) 
Cree a) 


where e; is the cross axis error term and Z is the error in the acceleration due to the 
cross axis sensitivity. 
2. Results and Validation 
Accelerations were measured for step aileron, elevator, and rudder inputs, 
and the results then compared to the actual accelerations computed for the equations 


of motion for the aircraft. A linear synthesis model was used for the initial testing, 
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while the results from the non-linear model were used for validation of the accelerom- 
eter model. Figures 5.5, 5.6, and 5.7 show comparisons between measured and actual 
accelerations generated from simulations of the non-linear model. The accelerations 
computed for the longitudinal and lateral cases were in close agreement with the 


accelerations computed by the non-linear model. 
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Figure 5.5: Measured and True Acceleration From a Step Elevator Input 
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Figure 5.6: Measured and True Acceleration From a Step Aileron Input 


Actual Thrust Acceleration, Y-Axis —— 
Measured Acceleration, Y-Axis °---~ 


Acceleration, g's 





Time, seconds 


10 


Figure 5.7: Measured and True Acceleration From a Step Rudder Input 
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B. RATE GYRO MODELING 

The rate gyros consist of a rotating disk, mounted in a gimbal mechanism. 
Though both single and two degree of freedom gyros are common, the tri-axial angle 
rate gyro supplied in the IMU-600D is modeled as three single-degree-of-freedom 
rate gyros, each measuring the angular rate along a particular axis. The dynamics of 


a gyro can be modeled [Ref. Sil 91] using Euler’s law 
EON a8 be 


wher {B} and {G} are the body and gyro coordinate systems, respectively. This can 
be expanded to 


d 
aN: = ee x sp Oe + GR—° Lo 


where °Lg = Ig?weg and the time derivative, £61 G is zero When the wheel rotates 
at a constant speed. Thus, except for the period when the wheel is coming up to 


speed, the equation for gyroscopic motion in a rate gyro can be written as 
aA = Bie x ip (Sao) 


where ? Ng is the torque vector acting on the gyro element and ?we is the angular 
rate of the gyro frame. Transfer functions can be developed for the torque input to 
the input axis as shown in Figure 5.8 and the rate output of the output axis. This 
derivation was not developed since data required for the gyro disk inertia, /g and the 
speed of rotation, wg, among other terms, is not available from the manufacturer. 
The rate gyros were modeled as second order transfer functions, with w, = 50 Hz. 
A third order Chebychev filter with a cut-off frequency of 20 Hz was added to the 
gyro model as an anti-aliasing filter. This Chebyshev filter was identical to the ones 
described for the accelerometers. The block diagram of the rate gyro model is shown 


in Figure 5.9. The linear synthesis model is shown in Figure 5.10 and was used to 
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Figure 5.8: Functional Diagram of a Rate Gyro [Ref. Bro 64] 


test the gyro models. In Figure 5.9 the gyro elements are shown for each axis and 
blocks for the error calculations are shown as well. 
1. Error Modeling 

Error terms are also present in the rate gyros and are due to either physical 
location on the aircraft, or mechanical errors, as in the accelerometers. Errors due 
to physical placement away from the c.g. can be corrected by using the equations 
derived for the linear and angular acceleration of an arbitrary point on the aircraft, 
shown later in this chapter. Mechanical errors are defined in the same manner as was 


done for accelerometers in Section A, and are listed in Table 5.2 along with other 
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Figure 5.9: Rate Gyro Model 


TABLE 5.2: GYRO CHARACTERISTICS 


[__Rotational Rate Range | £114. bdeg/sco_| 
[Rotational Rate Bandwidth |__20Hz | 
[Rotational Rate Bias__| 2.0% of Full Seale | 
[Rotational Rate Scale Factor | 0.5% of Full Seale | 
[ Rotational Rate Noise Floor | 0.05% of Fall Scale | 
[Gross Axis Sensitivity | 0.5% of Full Seale | 








important characteristics. The cross axis error is modeled in the same way as for 
the accelerometers by the use of the same matrix for computing errors in the angular 
rates. 
2. Results and Validation 
Angular rates were measured for step aileron, elevator, and rudder inputs 
and compared to the actual angular rates computed from the equations of motion. 


First a linearized synthesis model was used in the testing stages; then the sensors 
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Figure 5.10: Rate Gyro Synthesis Model 


were integrated with the non-linear aircraft simulation. Comparisons of actual and 
measured acceleration are given in Figures 5.11, 5.12, and 5.13. These figures show 
the accuracy of the angular rate sensors in measuring the angular rates computed 


with the non-linear equation of motion model. 
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Figure 5.11: Measured and True Angular Rates From a Step Aileron Input 
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Figure 5.12: Measured and True Angular Rates From a Step Elevator 
Input 
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Figure 5.13: Measured and True Angular Rates From a Step Rudder Input 
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C. PITCH, ROLL, AND HEADING SENSOR MODELING 

The last group of sensors to be modeled are those which measure the pitch, roll, 
and heading angles. The pitch and roll angle measurements are made with liquid 
pendulous devices. These are devices that have an electrolyte contained in a vial, 
and by measuring the capacitance changes of the vial as the electrolyte moves in 
response to aircraft angle changes, the pitch and roll angles can be determined. The 
heading sensor is a magnetic heading sensor. 

The primary concern with angular measurements is that accurate readings are 
obtained, even when the aircraft is experiencing linear and angular accelerations. 
When these errors cannot be corrected, the control system must be able to compen- 
sate for the errors. An inclinometer is typically modeled as a pendulum, shown in 
Figure 5.14 attached to a block, which can be considered to represent an aircraft. 
The equations describing the motion of the pendulum are derived in detail later in 
this section. The transfer function for the pendulum inclinometer can be represented 


as a second order transfer function with w, = 0.8 Hz and ¢ = 0.5. [Ref. WAT 93] 


O 5.03? 


© s?4 50.38 + 5.03?” 
1. Error Modeling 
The inclinometers are subject to several sources of errors, from both linear 
and angular accelerations, and from mechanical imperfections. Mechanical errors are 
listed in Table 5.3, and modeled as shown in Figure 5.15. Errors due to angular 
velocity can be compensated for with a complementary filter. Selecting the time con- 
stants appropriately will allow direct measurements of the angle from the inclinometer 
at low frequencies, and from integrated angular rates at higher frequencies where the 


inclinometer is inaccurate. 





TABLE 5.3: INCLINOMETER AND HEADING SENSOR CHARAC- 


TERISTICS 


It was determined that the effects of linear acceleration were also important 
to model. In order to model these effects it was necessary to derive a transfer function 
from the aircraft’s linear acceleration to the angle caused by that acceleration. The 


starting point was to model the inclinometer as shown in Figure 5.14. The derivation 


Figure 5.14: Simple Pendulum Inclinometer 
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Figure 5.15: Inclinometer Error Modeling 


results in two coupled equations, the first of which will be ignored, since the equations 
of motion for the aircraft are known and the motion of the inclinometer should have 
little effect on the aircraft motion. The second equation is for the motion of the 
pendulum as influenced by the aircraft’s acceleration and is what was used for the 
linear acceleration errors. 

First it was necessary to define the coordinates used to describe the incli- 
nometer system as g = (7,6) (see Figure 5.14). Next, since Lagrangian methods were 


used for the derivation, the kinetic energy, T, of the system was defined as 


T = 1/2Mz? +1/2m(z + 16cos6)? — 1/2m(16 sind)? 


= 1/2Ma’ | Woke i 2 (5.6) 


where P, and Py, represent the position of the pendulum in Cartesian coordinates. 


The potential energy, V, of the pendulum can be written as 
V =V,, = mgl(1 — cos8), (ou) 


where M, m, and ! are the mass of the aircraft, mass of the pendulum, and distance 
from the aircraft c.g. to the pendulous mass. The terms,P; and @, represent the 
position vector to the pendulous mass and the angle made by the pendulum with the 


vertical plane. The governing equation for Lagrangian dynamics is 
oper (5.8) 


where £ is the Lagrangian operator, £ = T — V, and Q, represents the non- 


conservative forces acting on the body. After some rearranging, 
L=1/2(M + m)z? + mlz6 cos + 1/2ml*6? — mgl(1 — cos8) (5.9) 


The result in Equation 5.9 can be substituted into Equation 5.8 resulting in 


aL 
Ox 
OL 
ag 
aL 
Ox 
aL 
BY, 
4 o£ 
dt Oz 
doe 
dt 06 


= (M+m)z + ml6cos6 

= mltcosé + 16 

= 0 

= —mlzé sind — mgl sind 

= (M+4+m)#+mlcos66— ml sind 
= mlcos6z — mlsin6 6 + ml’6. 


When these partial derivatives are substituted into Equation 5.8, the resulting equa- 


tions of motion for the pendulum are 


(M + m)é + mlcos0 6 — ml sind 6? 


I 
= 


Qo, 


ml cos6z + ml?6 + mglsin@ 
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where @, and @2 are terms representing damping in the system. The term Q, repre- 
sents damping on the aircraft for Q; = —(,z and the term Q»2 represents the viscous 


damping of the pendulous mass, Q2 = —3,0. The equations can be linearized, where 


| 
= 


(M+m)z+ml6+ 8,2 (5.10) 


I 
= 


mlé + m6 + 3,0 + mgl6 (5.11) 


Now it is apparent that Equation 5.10 is completely determined by the aircraft equa- 
tions of motion that have already been modeled. A Laplace transform of Equa- 


tion 5.11 can be performed to find the desired transfer function 


8) —l/l 
oe — soa (5. ap 
G gt eo Oye 





where the similarity to the standard second order transfer function is noted, making 
the term g/l = w?2. The term —1/] in the numerator can be solved for by substituting 
g = 32.174 f/s? and w, = 5.027 rad/sec, or 0.8Hz, with the result 1 = 1.27 jae 
result is then substituted into the block diagram for the inclinometer models shown 
in Figure 5.15. Responses for step inputs are shown in Figures 5.16, 5.17, and 5.18. 
In the lateral cases, there is very little difference between the measured and actual 
angles. The longitudinal case is quite different. In Figure 5.16, there is a considerable 
difference between the actual pitch angle, O, and the pitch angle measured with the 
inclinometer model. This difference the linear acceleration of the aircraft as a step 
elevator input is applied, and must be compensated for before a reliable control system 


can be developed. 
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Figure 5.16: Response to a Step Elevator Input 
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Figure 5.17: Response to a Step Rudder Input 
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Figure 5.18: Response to Step Aileron Input 
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The error in inclinometer readings due to linear acceleration has been de- 
termined experimentally by the manufacturer, and can be represented as a second 


order transfer function from g’s to O, 


O ; 157.087 


= "9 ) 1E7 NO. 1127 No2? Bh’ 
Pn ee mee 0s (5.13) 


where the gain, A’, was determined from an input of 1.2 V at 10 mV/g, and the 


resulting output of 0.4 V at 60 mV/deg. Thus the gain, A’, was 17.99. 


D. MODELING OF AN ARBITRARY SENSOR PLACEMENT 

The actual sensor placement in either the Bluebird or the AROD aircraft is not 
at the c.g., as was assumed in the previous sections. In order to model the actual linear 
and angular accelerations at the sensors, regardless of the position on the aircraft, 
new equations of motion must be derived. These equations can then be applied to 
the sensor inputs to obtain a proper output from the sensor. First the equations 
for linear acceleration will be derived using Newton’s third law for conservation of 
momentum, then the equations for angular acceleration will be derived using Euler's 
law for conservation of angular momentum. These equations must be expressed in 
terms of {B} since all the information measured by the sensors will be determined in 
the body coordinate system. 

1. Linear Accelerations 

In the derivation for equations of motion of an arbitrary point on the air- 

craft, the coordinate systems representing the inertial and body coordinate systems, 
similar to those shown in Figure 3.1 are used. Supplementing the derivation of equa- 
tions of motion for an aircraft, the motion of the sensor location on the aircraft, given 
by Po is necessary. In the inertial coordinate system, {U}, the position of Q can 


be written as “Po = 3R®Po9 +” Po. The velocity is first determined by applying 


(7 


Coriolis’ theorem from Equation 3.24 to obtain the first derivative. Here, 


"Vo = U Po = ” Pgo alg VS (E RP) (5.14) 
where 
d d 
5 (BRPa) = 85 (gRPo) + we x (BR*PQ). (5.15) 


The terms Ud and Bd refer to derivatives taken with respect to the inertial and body 


coordinate systems, respectively. The resulting expression for velocity is 
Vo = "Veo +'% x (gR”Po), (5.16) 


where ?£(2R®Pg) = 0 since Q is a point fixed on the aircraft. Accelerations are 


derived by applying Equation 3.24 twice to Equation 5.16 
ie Beda 
"Vo = ("Veo + & x (BRPQ)) +3 x ("Veo +B x (BR™Po)). (5-17) 


The BZ (-) term is differentiated, resulting in 


d d S 
2 ("Veo +% x (gR°PQ)) = ° = ("Veo) + Vup x RPA +3 x° Vo. (5.18) 


Substituting the results from Equation 5.18 into Equation 5.17, the result is 


d ie | | 
“Vo = °= (Veo) + 2 “we x UVo + Mbp x ZROPo +4 x (3 x GR Po) oe 


The desired result in {B} is obtained by simply premultiplying Equation 5.19 by BR 


dees 
°Ve = 5 (Puno) + 2 Pup x °VQ + Pup x PQ + Pwe x (Pup x “PQ), (5.20) 


where the identity [Ref. Sp 89] 
i R(wp x “ Vo) = (7 Rws) x (GR° Vo) 


was used. This result can now be substituted into Equation 3.27 and the resulting 


expression can be solved for £(? 


; UBo). 
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2. Angular Accelerations 
No further work is needed to derive the expressions for angular velocity 
and acceleration at a arbitrary point, since for a rigid body, these quantities remain 


the same anywhere on the body. 


WBO = WO VQ € { B} 


WBO sa WO VQ e {B}, 


Since w X w = 0. 


OS 


VI. CONCLUSIONS AND 
RECOMMENDATIONS 


A. CONCLUSIONS 


Based on the data presented in this thesis, the following conclusions are made: 


High-fidelity models of both the AROD and Bluebird aircraft were implemented 
using SIMULINK software. Use of the block diagram structure of the model 
allows changes to be easily made and requires no programming ability, other 


than the use of MATLAB. 


These models accurately represent the sensors, actuators, and engines associated 


with the particular aircraft. 


Expected errors in the accelerometers and rate gyros are very small. Pitch 


errors from the inclinometers, due to linear accelerations, will be much greater. 


The models that were developed only represent one flight condition. The AROD 
was modeled only in a hover and the Bluebird was only modeled in a cruise 
condition. Further work on the control, guidance, and navigation systems will 
be concentrated in these flight regimes. The data files for a given flight condition 


are easily replaced with tables, when more test data are available. 


B. RECOMMENDATIONS 


Based on the conclusions presented above, and the problems encountered during 


this study, the following recommendations are made. 


e The Department of Aeronautics and Astronautics should update the UNIX labs 


to include SIMULINK with MATLAB 4.0. This will allow increased instructional 
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use in flight dynamics, controls, and avionics courses. Additionally, a practical 


research tool would be available. 


The work in this thesis must be modified to include the aerodynamic charac- 
teristics of the Archytas aircraft. Achieving this step will require wind tunnel 
testing for the Archytas in all the expected phases of flight, especially the ex- 


tremely non-linear transition phase from vertical to horizontal flight. 


Integration of the quaternion based rotation matrix should be completed in 
order to take advantage of the superior qualities of quaternions in computational 


models. 
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APPENDIX A: MATHEMATICAL 
PROPERTIES 


In this appendix, some of the mathematical properties used in the text are 


described. 


A. CROSS PRODUCT PROPERTIES 
In this section, important properties of the cross product and cross product 


matrices are described The cross product between two vectors, 


ar b 
A= a, || "and Ss a\eoe 
b, 


a, 
a,0, er G30; | 





is defined by 


Ath. 





G70, — a0; 
Properties of the cross product are: 


0 —a, dy 
A <= {Ace wheter a, 0 —a, 


and is called a cross product matrix 


eo 2R(V xU)=(SRV) x (4RU) if ZR is a rotation matrix with V and U in the 
same coordinate system. That this matrix distributes across the cross product 


is obvious since rotation matrices preserve space geometry. 
e 2R(V -U) =(SRV) - (gRU) for the same reasons as above. 


oAx(Bx Cl] (45ers 


OP 
Le) 


Wie Cy bX (Ax C)+C x (Bx A) 
eax) ——b x A 
e -Ax =A’x 


B. DERIVATIVES OF VECTORS 
For any free vector, ?Vo (i.e. a velocity vector and any rotation matrix, #R, the 
derivative of the velocity of Q computed in {B} and expressed in {A} and denoted 


as “(7 Vo) is given as [Ref. Sil 91] 


d 
—(4(?Vo)) = = (Bh © Vp) 
ae a 
= BR Va + BR (PVG) 
: d 
Sy (in AR)BR EVE + BR(PVQ) 


since as shown in [Ref. Cra 86], 


d d 
<(4(7VQ)) = 4p x 4(FVQ) + AR= (VQ) 


so then 
APBp_A 
The same process can be carried out for an angular velocity vector 40g, 


5 (405) = 2(§R (408) 


dt 
d 
= “Op x “Op + pR—(°("Ox)) 
d 
= BR—(°(“08)) 


And if the origins of the coordinate systems {A} and { B} are coincident, the derivative 


can be expressed as 


For an applied vector ?Po (i.e. the position of point Q in the {B} coordinate system) 
and a rotation matrix #R, the time derivative of the position vector of Q, expressed 
in {A} , $( (4 Pg) as a function of it’s derivative computed in {B} is given by [Ref. 
Sil 91] 


d d 
= (“Pa) = = (BR @Po +“ Pao) 
d 
= BR MPa + BR (Pq) + “Ve 


= AQ, x (AR 2Po) +4R © (Po) +4Vp 
Therefore the velocity of the point Q can be expressed in {A} as 
AVq = “Op x (BR PQ) + BR OVO + “Ve 
or expressing the velocity of Q in the {B} coordinate system 
© (4Vq) = 7 (208) x °Pa + ?Vq + °(4VQ) 


In the case where the origins of {A} and {B} are coincident then the resulting ex- 
pression 1s 


Bug = Pwex 9Po + ?Vo + Fup. 


C. EQUATIONS USED FOR LINEARIZATION 
For any vector equation, H(c) = a x b where c = [a 5] and a, b, and c are all 


vectors, the Taylor series approximation can be written as 








OH 
it (@) = H (co) + Be lexeo Oe + eG 
and 
OH _ OH OH 
Ae © ne —|a= nyo = Ab leapged 
Now 0H/da can be written as 
OH  Oaxb)  O(-bxa)_ A-6) da 7 
da Ga Oa Oa X 0 — 0 ea ox 


Using the same relations, 0H/06 can be written 


OH _ O(-bxa) 
Ob Ob 


= ax 
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APPENDIX B: NUMERICAL RESULTS 


A. AROD RESULTS 


The following results are from the numerical linearization of the kinematic equa- 
tions of motion for the AROD. Using the trim command, based on an initial Op = 7/2, 


resulted in the state vector and input vector of: 


and u = 


a ey 


So 2 2 a & 


& 
| 
— 
qr 
~] 
=) 
S.> oo.o. SS OS  es 


The linmod command, [a,b,c,d]=linmod(’basInmod’,x,u), produced the following lin- 


earized system: 


OOF 6 0 0 
000 0 0 0 
pe 0 Ot 0 0 
O70 Ome 0 0 
000 0 0 —1.5024 
0 FOG) 10> or 0 | 
1 0 0650-70 
0 150 come 
,-{9 9100 0 
Coe ie CO 
00001 0 
b020 0 0 Oey 


and the c and d matrices were empty, since no outputs were defined. 
The following results are for the numerical linearization with gravitational ef- 


fects added to the 2-3-1 Euler angle kinematic equation model. The trim com- 
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mand, using the same initial conditions as before, resulted in the following vectors: 


[x,u]=trim(’grav]mod’,x0,|],|[],ix,[],[}) 


—0.00 
Dan 32.1740 
—0.00 
—0.00 
0.00 
—0.0256 
c= 0.00 anit. 
—(.00 
0.00 
0.00 
0.00 0.00 
1.5700 , 
0.00 


Linearizing the system, based on the state and input vector found by trimming at 


the desired conditions resulted in: [a,b,c,d]=linmod(’gravlmod’,x,u) 


0 0.00 0 0 0.00 0.00 0 -—0.0255 0.0002 
—0.00 0 0.00 —0.00 0 0.00 0.0256 —0.00 32.1740 
0.00 —0.00 0 —0.00 —0.00 0 -—0.00 —32.1740 0 
0 0 0 0 —0.00 —0.00 0 0 0 
i 0 0 0 0.00 0 —1.5204 0 0 0 
0 0 0 —0.00 1.5112 0 0 0 0 
0 0 0 1.00 —0.00 0.00 0 0 —0.00 
0 0 0 0 1.00 —0.00 —0.00 0 0.00 
0 0 0 0 0.00 1.00 —0.00 0 0 
and 
1.00 0 0 0 0 0 
1.00 0 0 0 0 
ie 0 1.00 0 0 0 
0 0 


0 
0 
0 0 0 1.00 
0 
0 


0 0 OO O 1.00) 


Trimming the full EOM model at hover conditions resulted in the following: 


—(0.00 
—0.00 
0.00 
—0.00 
fo 0.00 | andu= 
—0).00 
0.00 
1.5708 
0.00 


0.00 
0.00 
0.1807 
6387.2 
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The subsequent linearization for the complete model vielded: 


0) 0 0 0 0 0.0002 
0.00 0 0 0.00 0 0.00 0.00 —0.00 
070 0 S000 0 0 -—0.00 —32.1740 
00 0 0 0 0 0 0 
= 0 0 0 —0.00 0 -—1.5174 0 0 
0 0 0 -—0.00 1.5082 0 0 0 
00 0 £1.00 —0.00 0 —0.00 0 
0 0 0 0 1.00 —0.00 0.00 0 
0 0 0 0 0 1.00 0.00 0 
and 
0 0 QO 0.0112 
0 0 0 0 
0 0 0 0 
0 Q 24.8194 0.0003 
b= | —6.6192 0 0 —0.00 
0 -—6.5791 0 —0.00 
0 0 0 0 
0 0 0 0 
0 0 0 0 
Computation of the eigenvalues gave: 
0 
0 
0 
0 
eig(a) = 0 
—0.00 
0+ 1.51282 


0 — 1.51287 


0 


0.0002 
32.1740 


Now consider the quaternion based model Rpm is fixed at 6800 RPM. 


To hold the desired vertical attitude, the quaternion states were held fixed: 


8S 


ee SS) ee SS 


where the state vector of nominal conditions was: 


Le 


Trim and linearization of the quaternion-based gravitational model resulted in the 


following state and input vectors, as well as linearized a and b matrices: 


—0.00 
0.00 
—0.00 32.1740 
0.00 —0.00 

~ _ | 70-00 | 0:00 
0.00 | ° 0.00 } ° 
0.7071 0.00 
—0.00 0.00 | 
0.7071 
0.00 | 

a =, columns 1-6, 

0 0.00 0.00 0 0.00 0.00 
—0.00 0 0.00  —0.00 0 0.00 
—0.00 -0.00 0 -—0.00 —-0.00 0 

0 0 oO 0 0.00 —0.00 

0 ee a0 0.00 0 —1.6051 

0 OF - 10 0.00 1.6062 0 

0 0 0 0.00 —0.3536 —0.00 

0 (em 0 3536) —0.00) 0.3536 

0 0 0 0.00 0.3536 0.00 

0 0 O 0.3536  —0.00 0.3536 | 
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columns 7-10 


—45.5009 0 —45.5009 0 
0.00 —45.5009 0.00 45.5009 
45.9012 —0.0003 -—45.5012 0.0003 
0 0 0 0 
0 0 0 0 
0 0 0 Ce 
0 —0.00 0.00  —0.00 
0.00 0 0.00 0.00 
—0.00 —0.00 0 0.00 
0.00 —0.00 —0.00 0 
and 
1.00 0 0 0 0 0 
0 1.00 0 0 0 0 
= 0 0 1.00 0 0 0 
0 0 0 1.00 0 0 
0 0 0 0 1.00 0 
Pe 0 0 0 0 1.00 | 
Eigenvalues are computed as: 
—0.00 + 0.002 
—0.00 — 0.002 
0.00 
0 
eae —0.00 + 1.60572 
—0.00 — 1.60572 
—0.00 
0.00 
—0.00 + 0.002 
—0.00 — 0.002 


B. BLUEBIRD NUMERICAL RESULTS 


Results for the numerical trim and linearization of the Bluebird model are given 


below. The results are from the trim of the kinematic equations of motion for the 


90 


Bluebird model are based on the nominal conditions of: 


72.9954 
0 
6.6757 
0 


© 


r= 


Oo © 


0.0912 
0 


with the states, u, w, and O held constant, zr = [1 5 8]’. The trim command, 
(x,u}=trim(’basic’,x0,|],|],1x,[],[]) resulted in the following state and input vectors: 


12.9954 
0 
6.6757 
v 


5 Guile, Ye — 


=) 
II 
© 


(me i = 
2S] 2 2 2 © 


NEOOT2 LO 
0 


The numerical linearization of the basic model, [A,B,C,D]=linmod{(’basic’,x,u), re- 


sulted in the following A and B matrices: 


00 0 De 6.8750 0 
0 0 0 6.6757 0 —72.9954 
pa | oot 0 72.9954 0 
00 0 0 0 0 
00 0 0 0 0 
000 0 0 0 | 
and 
Io? 0 0 O @ 
De Oe 
Oe Ona0n.G 
= li OO i 0 a 
000010 


000001] 


The C and D matrices were empty since no outputs were defined. 
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Results for trim and linearization of the kinematic model with gravitational 


effects added yielded the following state and input vectors: 


72.9954 
0.00 se 
6.6763 ec 
0.00 0.00 
—31.4605 
i 0.00 |, andu= 
9.00 —0.00 
0.3945 
—0.00 | ~0.00 | 
0.0912 
0.00 


The numerical linearization of the model resulted in the following A and B matrices: 


A= 
—0.0007 -—0.00 0.0079 0 -—6.6763 0.00 —0.00 -—32.0451 0 
0.00 0 0.00 6.6763 0 —72.9954 32.0403 0.00 0 
—0.0001 0 0.0007 -—0.00 72.9954 0 -—0.0002 —-—2.8773 0 
—0.00 0.0087  —0.00 0 —0.00 0.00  —0.00 —0.00 0 
—0.00 0.00 0.0005 0.00 0 —0.00 0.00 0.0361 0 
—0.00 0.0010 0.00 —0.00 —0.00 0 —0.00 —0.00 0 
0 0 0 1.00 —0.00 0.0915 0.00 —0.00 0 
0 0 0 0 1.00 0.00 0.00 0 0 
0 0 0 0 —0.00 1.0042 0.00 —0.00 0 
and 
L000 Foe 
0100 0 0 
0 0 105080 
is 00010 0 
OE On Oe Gn 
:0 0 0 0 05 J 
Results for trimming the entire model were A = 
—0.0622 0.00 0.3431 0 —1.6187 0.00 0.00 -—32.0416 0 
0.00 —0.3865 0.00 1.7450 0 —71.6817 32.0403 0.00 0 
—0.7558 0.00 —4.7145 0.00 67.1233 0 -—0.0002 —-—2.8771 0 
0.00 —0.1455 0.00 —5.3700 0.00 1.5003 0.00 0.00 0 
0.0155 0.00 —0.1907 0.00 —3.1266 0.00 0.00 0.0362 0 
0.00 0.1418 0.00 —1.0589 0.00 —0.7970 0.00 0.00 0 
0 0 0 1.00 0.00 0.0915 0.00 0.00 0 
0 0 0 0 1.00 0.00 0.00 0 0 
0 0 0 0 0.00 1.0042 0.00 0 0 
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and 


—4.3850 0 Q 8.7745 

0.00 5.6803 0 0 

—37/.8929 0 0 0 

0.00 0.6216 45.9851 0 

a2 eo 0 0.00 0 
0.00 —7.1281 —6.1465 0 

0 0 0 0 

0 0 0 v 

0 0 0 0 


Complete model linearized with Og = 0 instead of Og = ao The initial conditions 


are Now: 


73.3000 
0 
0 
0 
C0 0 
0 
0 
0 
0 
The state and input vectors obtained from trimming at this state are 
73.3000 
—0.00 
eee ~0.0181 
—0.00 
ie —0.00 | andu= Su 
0.00 0.00 
0.00 0.2336 
—0.00 
—().00 
with the linearized A and B matrices: A = 
—0.0635 0.00 0.3277 0 —1.4922 0 —0.00 —32.1740 
—0.00 —0.3911 —0.00 1.6086 0 -—72.6109 32.1740 —0.00 
—0.7572 —0.00 —4.7741 0 67.9934 0 -—0.0002 -—0.0002 
0.00 —0.1471 —0.00 —5.4414 —0.00 Ipols3 0.00 —0.00 
0.0151 —0.00 —0.1933 0 —3.1672 0 0.00 —0.00 
0.00 0.1440 0.00 —1.0578 0.00 —0.8114 0.00 —0.00 
0 0 0 1.00 0 —0.00 0 —0.00 
0 0 0 0 1.00 —0.00 —0.00 0 
0 0 0 0 0.00 1.00 —0.00 —0.00 
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Se SS OS Qe a a 


— 4.5835 0 0 8.7745 
0.00 5.8282 0 0 
—38.8681 0 0 0 
—0.00 0.6252 47.1717 0 
B= |} —22.0417 0 0 0 
—0.00 -—7.3151 —6.4345 0 
0 0 0 0 
0 0 0 0 
0 0 0 0 
The eigenvalues for the case where Og = 0 are 
0 
—0.5285 + 3.63462 
—0.5285 — 3.63462 
—5.6291 
erg(A} = | —3.9833 + 3.55212 
—3.9833 = dido2h 
0.0420 


—0.0191 + 0.49632 
—0.0191 — 0.49632 


C. CESSNA 172 RESULTS 


The Cessnal 172 model was trimmed and linearized using the state vector 


Se 
pt 
ce) 


Bay 
= 
| 
SAS SS. O Sre 


as specified in [Ref. Ros 79]. The states u, w, and O were held constant making the 
term go 193 312 


Only the complete model was linearized, with the results of the trim routine 


Q4 


given as: 


—0.0442 
0.00 

mall co10 
0.00 
0.0024 
—0.00 

0 

0 

0 


and 


and the following eigenvalues 


219.00 
0.00 
he 0.0001 
—0.00 
je —0.00 | andu= HO 
0.00 
See 1.0013 
—0.00 , 
—0.00 
0.00 
The linearized model had the following results for the A and B matrices A = 
0.00 0.0848 0 0.0382 0 
—0.1620 —0.00  —0.0394 Q —217.2141 
—0.00 —2.1805 0 212.5399 0 
—0.1313 0.00 —12.4093 0.00 Peale 
—0.00 —0.1085 0 —6.0778 0 
0.0462 0.00 —0.3807 0.00 — 1.2600 
0 0 1.00 0 —0.00 
0 0 0 1.00 —0.00 
0 0 0 0.00 1.00 
—6.2509 0 053.2259 
—0.00 19.4571 0 0 
—44.3392 0 0 0 
0.00 4.7446 57.4954 0 
B= | —39.4919 0 0 0 
—0.00 —10.2288 —8.2562 0 
0 0 0 0 
0 0 0 0 
0 0 0 0 
0 
— 12.4309 
—0.6947 + 3.30802 
—0.6947 — 3.30802 
e7g(A) = | —4.1303 + 4.38952 


—4.1303 — 4.38952 

—0.0109 
—0.0209 + 0.17942 
—0.0209 — 0.17942 
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0.00 
32.1740 
—0.0002 
—0.00 
0.00 
0.00 

0 

0.00 
—0.00 


—32.1740 
—(0.00 
—0.0002 
—0.00 
—0.00 
—0.00 
—0.00 

0 

—0.00 


(em) (as as | as oe) es) ae ae) 


APPENDIX C: PROGRAM LISTINGS 


A. AROD MATLAB ROUTINES 
1. Main Routine 


function accel=main(vstate,m,rho,A,R) 


4 Function will compute the accelerations due to the 

h gravitational forces, aerodynamic forces and moments, 
h and control forces and moments. 

Zi: The values for S,rho,m,b,and c are used as input 
if arguments to the function call, and are loaded 

y/ from the workspace. There should be a file, 

h loaddata.m loaded prior to running the 

h simulation. 


7, define states in terms of the input vector 
u=vstate(1); 

v=vstate(2); 

w=vstate(3); 

p=vstate(A4) ; 

q=vstate(5) ; 

r=vstate(6) ; 

phi=vstate(7); 

theta=vstate(S8); 


psi=vstate(9) ; 
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4, Define the control inputs 
de=vstate(10) ; 
dr=vstate(11); 
da=vstate(12); 


drpm=vstate(13) ; 


% this subroutine computes linear and angular accelerations 

% given angular and linear velocities; 

% the input is 6x1 vector = [uvwpdqr]’ 

%, the output is feedback part of d/dt [uv wpqr]’ 

% the output also includes the Euler angle derivatives, based on 


% a 2-3-1 transformation, for Ldot, used in the function grav.n. 


v = vstate(1:3); 


omega = vstate(4:6) ; 


vdot = -crpr(omega,v) ; 

tib,ir] = inertia; 

% compute the angular momentum due to the body’s inertia, Ib 

Lb = Ib *omega; 

4 compute the angular momentum due the spinning prop’s inertia, Ir 
OmegaR=drpm*2*pi/60; % angular velocity of the prop, rad/sec 
Lr=Ir*[OmegaR;0;0]; 

temp=Lr+Lb; 


omdot = -Ib\crpr(omega, temp) ; 


oT 


vstatedot=[vdot;omdot] ; 


4 Use the Euler angle propagation equations for a 2-3-1 rotation sequence 
Ldot=[p/cos(psi)-cos(phi)*sin(psi) /cos(psi)*qtsin(phi)*sin(psi)/cos(psi) *r; 
cos(phi)/cos(psi) *q-sin(phi)/cos(psi)*r; 


sin(phi) *q+cos (phi) *r] ; 


h Given the vector containing the state derivatives, 

ih The function will compute the forces and moments due to 
ts the control derivatives, Cfd, where this 

h is det/dd. 

Wf The values for rho,A,R,m are used as input 

h arguments to the function call, and are loaded 

h from the workspace. There should be a file, 

h arod.mat loaded prior to running the 

| simulation. 


7% hover case V=0, dimensionalize the control derivatives based on 
7 induced velocity through the rotor disk, Vi=sqrt(T/2/A/rho) 
%4 Calculate the quantities needed for the force 


% pealeulation: 


% Call a function to return the stability derivatives wrt to moments 
% Could put a call to a lookup table here 


74 Syntax--Z = TABLE2(TAB,X0,Y0) or Y = TABLE1(TAB,X0) 


Cfd=getcfd; 

%, Define the derivatives 

Clda=Cfd(4,3); 

Cmde=Cfd(5,1); 

Cndr=Cfd(6,2); 

% Calculate the Force due to Cfd derivatives 
74 No Aerodynamic Forces in a Hover 

D=0; 

Y=0; 


L=0; 


7% calculate the force due to thrust in body coordinates 


% USING THRUST VS. RPM FROM BOB STONEY’S TEST RUNS 


T=0 .0297*drpm-104.7; 

Vi=sqrt (T/2/A/rho) 

qbar=.5*rho*Vi"2;  % Vi is induced velocity, not forward speed 
Fout=[D;Y;L]; 


Fout=(Fout+[T;0;0])/m; 


% Calculate the Moment due to Cfd derivatives 

% ltr relates the duct swirl to the moment, 1, produced by thrust 
itr=-0 .0542*T-0 .9138; 

l=qbar*A*R* (Clda*da)+ltr; 

m=qbar*A*Rx* (Cmdexde) ; 


n=qbar*A*R* (Cndr*dr) ; 


oh) 


Nout=Ib\{l;m;n]; 


FNcfx=[Fout ;Nout] ; 


y/ Given the vector containing the state derivatives, 
h and euler angles, the function will compute 
h the forces due to gravity acting on the aircraft. 


y/ 

h 

% Calculate gravitational force, based on a 2-3-1 Euler angle 

4% rotation for position of the aircraft. Rotation matrix is Ru2body 


4 Z for {U} is positive down and X for {B} is straight up. 


FNgrav=32.174*[-sin(theta)*cos(psi) ; 
sin(theta)*sin(psi)*cos(phi)+sin(phi)*cos (theta) ; 
cos(phi) *cos(theta)-sin(theta)*sin(psi)*sin(phi) ; 


0;0;0]; 


% Sum up the normalized forces and moments to feed back into the integrator 


vstatedot=[vstatedot+FNcfx+FNgrav;Ldot] ; 


accel=vstatedot; 


2. Supporting Subroutines 


function y = crpr(omega,x) 


100 


=~ 


4 y = omega X x 


p = omega(1) ; 
q = omega(2) ; 
r = omega(3) ; 
a= (0 -r q; 
0 -p; 
-q p Ol; 
wat FX; 


function Fgrav=grav(x) 
% GRAV will compute the gravitational force 


% acting on the body, in body coordinates 


g=x(1); 
phi=x(2) ; 
theta=x (3) ; 
psi=x(4) ; 


Fgrav=g*[-sin(theta)*cos(psi) ; 


sin(theta)*sin(psi)*cos(phi)+sin(phi)*cos(theta) ; 


cos (phi) *cos(theta)-sin(theta)*sin(psi)*sin(phi) ] 


funetion Lib,ir] = inertia 


4, this subroutine creates inertia matrices 


101 


this subroutine computes the crossproduct of omega and x: 


called ib 


3 


7% and ir for the body and rotor inertia, respectively. 
% Arod hover case 
slip > ete De ZS hPa 


iyy = 3.9584; 


ZZ, S15 .9o7o 

uxz = 0; 

1rx=.00898; 

iry=.0045; 

1rz=.0045; 

ab = [ixx 0 -1xz:0 iyy)0;-i1xz 0 azz 


Pex 0 Oy 0 0 0 Oba ezir 


sD & 


3. Data and Initialization Subroutines 


¥, load bluebird data 


A=3.14- /, khrea of rotor dist ttac 
Vt=712.09; 7, Rotor tip speed, rad/sec 
m=2.6419; 4 Mass, slugs 

R=1; 74, Radius of Rotor Blade, ft 
rho=.002377; % Air density, slug/ft73 
Uo=0; 74, Nominal Velocity, ft/sec 


% Initial Euler Angle Orientation, radians 

Phi=0; 

Theta=1.57; 74 Same as Steady State alpha 
Psi=0; 


Lo=(Phi;Theta;Psil]; 


% Initial Conditions to determine the Aircraft 
% Linear and Rotational Velocity States 

% u,v,w are computed from Uo, alpha, and beta 
4% p,q,r are assumed as zero. 

alpha=Theta; 

beta=0; 


Xo=[Uo;alpha;beta] ; 


¥, Returns Initial Conditions for the main 
4 integrator in the rigid body EOM block 


i0=init_var(Lo,Xo); 


function Cfd=getcfd 

4 Cfd=getcfd will return values for 

4 The stability derivatives for D,Y,L,1,m,and n 
4 due to the control inputs. 

47, format for data is; 

% [CDde CDdr CDda 

peecide ... 

feecide ... 

% Clde Cldr Clda 

4 Cmde . 

famcnde ... Cnda 

%4 Data is non-dimensionalized by using induced velocity 


% Vissqrt(T/2/A/rho) 
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h 


for arod hover case 
Cfd=[ 0 0 0; 
0 0 Oo: 
0 0 0 
0 O 1.438; 
=12 25500 Cc. 


O -1.233 0]; 


function i0=init_var(Lo,Xo) 


h 


INIT_VAR(X) will initialize the integrators 
initial states, 10, given the initial 
conditions desired. 

Required initial conditions are the Euler 
angle orientation, total velocity, Uo, initial 
AOA, and sideslip angle, beta. 
Lo=[phi;theta;psi]’ 

Xo=[Uo;alpha; beta] ’ 


All body rotation rates are assumed to be zero 


Initialize the states u,v,w,p,q,r 


Vo=Xo(1); 


alpha=Xo(2) ; 


beta=Xo(3): 
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Ca=cos (alpha) ; 


Sa=sin(alpha) ; 


CB=cos(beta) ; 


SB=sin (beta) ; 


Rwb=[Ca*CB -Ca*SB -Sa;SB CB 0;Sa*CB -Sa*SB Ca]; 


v0=Rwb*[U0;0;0]; 


w0=(0;0;0]; 


70=([v0:wO:Lol; 


function Qo=initgq (lambda) 


h 
h 
y/ 
y/ 
h 


Function initQ will return values for 
the quaternion DCM based on a given 
set of Euler angles. 

set for a Euler 3-2-1 rotation 
Q(1)=BO 

Q(2)=B1 

Q(3)=B2 


Q(4)=B3 


phi=lambda(1); 


theta=lambda(2) ; 


psi=lambda(3) ; 


Qo(1)=.5¥*sqrt(1+cos(psi)*cos(theta)+sin(psi)*sin(theta)*sin(phi)... 
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+cos(psi)*cos(phi)+cos(theta)*cos(phi)) ; 


Qo(2)=1/4/Qo(1)* (cos (theta) *sin(phi)-sin(psi)*sin(theta)*cos(phi)... 


+cos(psi)*sin(phi)) ; 
Qo(3)=1/4/Qo(1)*(cos(psi)*sin(theta)*cos(phi)+sin(psi)*sin(phi)... 


+sin(theta)); 


Qo(4)=1/4/Qo(1)* (sin(psi) *cos(theta)-cos(psi)*sin(theta)*sin(phi)... 


+sin(psi)*cos(phi)); 
Qo=Qo’ ; 


B. BLUEBIRD MaTLAB ROUTINES 
1. Main Routine 


function accel=main(vstate,rho,b,c,S,m,Xo) 


h Function will compute the accelerations due to the 

Hf gravitational forces, aerodynamic forces and moments, 
y/ and control forces and moments. 

4, The values for S,rho,m,b,and c are used as input 
y/ arguments to the function call, and are loaded 

yA from the workspace. There should be a file, 

h loaddata.m loaded prior to running the 

h simulation. 


74, define states in terms of the input vector 
u=vstate(1); 
v=vstate(2) ; 


wevstate(3); 
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p=vstate(A) ; 
g=vstate(5) ; 
r=vstate(6) ; 
phi=vstate(7) ; 
theta=vstate(8) ; 


psi=vstate(9) ; 


4, Define the control inputs 
de=vstate(i0); 
dr=vstate(ii); 
da=vstate(1i2); 


dt=vstate(1i3); 


% calculate the aerodynamic terms 


Vt=sqrt (u*2+v"2+w 2) ; 


qbar=.5*rhox(Vt) “2; 


Ib=inertia; 


% wind to body transformation 


Rwb=rw2b(vstate, Vt) ; 


%4 CHI will compute the left hand side of the state 

4 equation. This is the term dependant on dCf/dxdot. 

4, Now calculate the S$ matrix that non-dimensionalizes 

4 the moments. Also includes the correction for Lift and Drag 
4 acting in the direction opposite to the positive coordinate 
4 direction. 

4 Get the stability derivatives for forces and moments 
Cixdot=getcfxd; 

CLad=Cfxdot (3) ; 


Cmad=Cfxdot (5); 


4 Twb is a intermediate step 


Twb=[eye(3)/m zeros(3) ;zeros(3) inv(Ib)]*([Rwb zeros(3) ;zeros(3) Rwb]; 


4, calculate the M2 matrix to allow use of the states, rather than 
4 the normalized states. Only accounts for the alpha dot variables 


4 since the beta dot terms are not ordinarily computed. 


4M2=[0 0 1 0 0 0] *c/2/Vt*2; 
ASprime=diag([-S S -S S*b S*c S*b]); 
4 To save some math here, the product of Sprime, Cfxdot, and M2 is: 
PROD=[0 0000 0;0 000 0 0;0 O -S*CLad*c/2/Vt"2 0 0 0; 
00000 0;0 0 Sxc*Cmad*c/2/Vt 2 0 070;0 0 Onc noRcle 
% Calculate chi 
chi1l=(eye(6)-Twb*qbar*PROD) ; 


h 


y | Given the vector containing the state derivatives, 
Z and euler angles, the function will compute 


h the forces due to gravity acting on the aircraft. 


% Calculate gravitational force, based on a 3-2-1 Euler angle 


4% rotation for position of the aircraft. Rotation matrix is Ru2body 
Fgrav=32.174*[-sin(theta) ; 
cos(theta)*sin(phi) ; 


cos(theta)*cos(phi)]; 


% Premultiply by the Chi*-1 term from the left hand side 


FNgrav=chii\[Fgrav;0;0;0]; 


4 Cfx(u) Given the vector containing the state derivatives, 


y/ The function will compute the forces and moments due to 
yA the stability derivatives, Cfx’, where this 
h is dCf/dx. 


h 
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% Call a function to return the stability derivatives wrt to moments 
4 Could put a call to a lookup table here 


% Syntax--Z = TABLE2(TAB,X0,Y0) or Y = TABLE1(TAB, XO) 


Cfix=getcfx; 

CDu=Cfx(1,1); CDa=Cfx(1,3); 

CYb=Cf£x(2,2): CYp=Ctx(2,4);>" CYr=Cix(Z.6). 
CLu=Cfx(3,1); CLa=Cix(3,3); CLlq=enx@ 7s): 
Clb=Cfx(4,2); Clp=Cfx(4,4); Clr=Cfx(4,6); 
Cmu=Cfx(5,1); Cma=Cfx(5,3); Cmg=Cfx(5,5); 


Cnb=Cf£x(6,2);  Cnp=Cix(6,4);  Cnr=cix(6. 6). 


SSqcerci 0; 
CDO=ss(1); 
CLO=ss(3); 


Cm0=ss(5); 


Cid=getcfd; 

4, Define the derivatives 
GNde=Cfd(1. 4). 

CYdr=Cfd(2,2); CYda=Cfd(2,3); 
CLde=Cfd(3, 1); 

Clidr=Cfd(4,2); Clda-erd(4na). 
Cmde=Cfd(5,1); 


Cndr=Cfd(6,2); Cnda=Cfd(6,3); 
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Pecatculate the Force due to Cfx’ derivatives 
% And the control derivatives 


% ain the wind coordinates 


D=-S*qbar/m* (CD0+CDa*w/Vt+CDdexde) ; 
Y=S*qbar/m* (CYb*v/Vt+CYr*r*b/2/Vt+CYdr*dr+CYda*da) ; 


L=-S*qbar/m* (CLO+CLa*w/Vt+CLq*q*c/2/Vt+CLdexde) ; 


% calculate the force due to thrust in body coordinates 


% THRUST IS ESTIMATED, BASED ON 4.0 HP, PROP EFFICENCY=.65 


T=15; 


Xt=T/m*dt ; 


Fout=Rwb*([D;Y;L]); 


Fout=Fout+[Xt;0;0]; 


% Calculate the Moment due to Cfx’ derivatives 

% And the control derivatives 

l=qbar*S*b* (Clb/Vt*v+Clp*b/2/Vt*p+Clr*b/2/Vt*r+Cldr*dr+Clda*da) ; 
m=qbar*S*c* (Cm0+Cma/Vt*w+Cmq/Vt*c/2*q+Cmde*de) ; 


n=qbar*S*b* (Cnb/Vt*v+Cnp*b/2/Vt*p+Cnr*b/2/Vt*r+Cndr*dr+Cnda*da) ; 


Nout=Ib\(Rwb*[{1;m;n]); 
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4 Premultiply by the Chi”-1 term from the left hand side 
FNcfx=chi1\[Fout ;Nout] ; 

h 

f wane nnn ene nnn ne 


4 this subroutine computes linear and angular accelerations 
I 
Y 
h 
i 


given angular and linear velocities; 


the input is 6x1 vector = [uvwpqr]’ 


the output is feedback part of d/dt [uvwpqr]’ 


the output also includes the Euler angle derivatives 


4 Ldot, used in the function grav.m. 


v = vstate(1:3); 


omega = vstate(4:6) ; 
vdot = -crpr(omega,v) ; 
temp = Ib*omega; 


omdot = -Ib\crpr(omega, temp) ; 


vstatedot=chil\[vdot;omdot] ; 


% Use the Euler angle propagation equations 
Ldot=[p+sin (phi) *tan (theta) *qtcos (phi) *tan(theta)*r; 
cos(phi)*q-sin(phi)*r; 


sin(phi)/cos(theta) *q+cos(phi)/cos(theta) *r] ; 


Ike 


vstatedot=[vstatedot+FNcfx+FNgrav;Ldot] ; 


accel=vstatedot; 


2. Supporting Subroutines 
function y = crpr(omega,x) 
% this subroutine computes the crossproduct of omega and x: 


4 y = omega X x 


p = omega(1); 
q = omega(2); 
r = omega(3) ; 
te—)10 -r q; 
me =p; 
=6| jo BIR 
yout * x; 


function Fgrav=grav(x) 
% GRAV will compute the gravitational force 


% acting on the body, in body coordinates 


g=x(1); 
phi=x (2) ; 


theta=x (3) ; 
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Fgrav=g*(-sin(theta) ; 
cos(theta)*sin(phi) ; 


cos(theta)*cos(phi)]; 


function Ib = inertia 

% this subroutine creates inertia matrix called Ib 
4 for the Bluebird test aircraft. 

% All units are slug-ft*2 

Ix=10.0; 

Ty=16.12; 

Iz=7 .97; 


Tb=14x 0.20; 0'-ly 07070 Iz: 


function [Rwb]=Rw2b(x,Vt) 


% RWIND2BODY Rotation matrix for wind to body coordinate transformations. 
Gaal Vt 
Sa=x(3)/Vt; 
CB=x(1)/Vt; 


SB=x(2)/Vt; 


Rwb=[Ca*CB -Ca*SB -Sa;SB CB 0;Sa*CB -Sa*SB Ca]; 
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3. Data and Initialization Subroutines 


% load bluebird data 


S=22.38; — % Planform Area, ft~2 
Uo=73.3; % Nominal Velocity, ft/sec 
m=1.7095; 4 Mass, slugs 

b=12.42; je opan, ft 

c=1.802; % Average Chord, ft 
rho=.002377 ; % Air density, slug/ft*3 


% Initial Euler Angle Orientation, radians 

Phi=0; 

Theta=.0912; % Same as Steady State alpha 
Psi=0; 


Lo=[Phi;Theta;Psi] ; 


% Initial Conditions to determine the Aircraft 
% Linear and Rotational Velocity States 

% u,Vv,w are computed from Uo, alpha, and beta 
% p,q,r are assumed as zero. 

alpha=Theta; 

beta=0; 


Xo=[Uo;alpha; beta] ; 


% Returns Initial Conditions for the main 
% integrator in the rigid body EOM block 


10=init_var(Lo,Xo) 
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function i0=init_var(Lo,Xo) 

% INIT_VAR(X) will initialize the integrators 

4 initial states, i0, given the initial 

% conditions desired. 

% Required initial conditions are the Euler 

4% angle orientation, total velocity, Uo, initial 
% AOA, and sideslip angle, beta. 

4 Lo=([phi;theta;psi]’ 

4 Xo=([Uo;alpha;beta]’ 


% All body rotation rates are assumed to be zero 


% Initialize the Euler angle DCM 


4 Initialize the states u,v,w,p,q,r 
Vo=Xo(1) ; 
alpha=Xo(2) ; 


beta=Xo(3); 


Ca=cos(alpha) ; 
Sa=sin(alpha) ; 
CB=cos (beta) ; 


SB=sin(beta) ; 


Rwb=[Ca*CB -Ca*SB -Sa;SB CB 0;Sa*CB -Sa*SB Cal]; 


vO=Rwb*[U0;0;0]; 
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w0=[0;0;0]; 


i0=[v0;w0;Lo]; 


function Cf0=getcf0 
%4 CfO=getcfO will return values for 
% the nominal values for coefficients 


%4 format of input is [CDO CYO CLO C10 CmO CnO]’; 


cf0=[0.03 0 0.300 0]’; 


function Cfd=getcfd 

%4, Cfd=getCfd_F(n) will return values for 
%4 The stability derivatives for D,Y,and L 
% due to the control inputs. 

4, format for data is; 

% [CDde CDdr CDda 

fecide ... 

% CLde .. 

% Clde Cldr Clda 

%4 Cmde ... 


yeecnde ... Cnda] 


% For the test aircraft Bluebird 


% Derivatives from DATCOM 


Lee 


Cfd=[ .065 0 0; 
0 .0697 0; 

472 0 0 
0 .0028 .265; 

-1.41 0 0; 


0 -.0329 -.0347]; 


function Cfx=getcfx 

4 Cfx=getcfx(n) will return values for 

4, The stability derivatives for D,Y,and L 
4, due to the state vector. 

7, format of data: 

((eDu-CDb CDba-Cbp eda Cir. 

y MCh Gh eens 

y eC) Gt Ota 

7 Ae Oa Roe 

i Griieean 


y alee Sugeet ake 1 


% For the test aircraft Bluebird 


Derivatives from DATCOM 


nd 


Cfx=[0 0 . 188 0 0 OF 
Omy0 31 0 0 0 JO0973= 


0 0 4.22 O 3.94 0; 


Ca 1co9i es 0 =]:2363 0 ale 
0 OPeat oS Oe La ( ee 0: 


Guen02s7 9 90) —.048) 20 - 0452]; 


function Cfxdot=getcfxd 

% Cfxdot=getcfxd(n) will return values for 
% The stability derivatives for D,Y,and L 
%Z due to the state vector. Beta dot terms are ignored 
% since they are not normally determined. 
% format is: 

%{ CDadot 

% CYad 

% CLad 

4 Clad 

74 Cmad 


paecnad | 


Cfxdot=[(0;0;1.32;0;-4.7;0]; 


NN: 
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